{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"# Tutorial: Gaussian pulse initial data for a massless scalar field in spherical-like coordinates\n",
"\n",
"## Authors: Leonardo Werneck and Zach Etienne\n",
"\n",
"# This tutorial notebook explains how to obtain time-symmetric initial data for the problem of gravitational collapse of a massless scalar field. We will be following the approaches of [Akbarian & Choptuik (2015)](https://arxiv.org/pdf/1508.01614.pdf) and [Baumgarte (2018)](https://arxiv.org/pdf/1807.10342.pdf).\n",
"\n",
"**Notebook Status**: Validated \n",
"\n",
"**Validation Notes**: The initial data generated by the NRPy+ module corresponding to this tutorial notebook are used shown to satisfy Einstein's equations as expected [in this tutorial notebook](Tutorial-Start_to_Finish-BSSNCurvilinear-Setting_up_ScalarField_initial_data.ipynb).\n",
"\n",
"## Python module which performs the procedure described in this tutorial: [ScalarField/ScalarField_InitialData.py](../edit/ScalarField/ScalarField_InitialData.py)\n",
"\n",
"## References\n",
"\n",
"* [Akbarian & Choptuik (2015)](https://arxiv.org/pdf/1508.01614.pdf) (Useful to understand the theoretical framework)\n",
"* [Baumgarte (2018)](https://arxiv.org/pdf/1807.10342.pdf) (Useful to understand the theoretical framework)\n",
"* [Baumgarte & Shapiro's Numerical Relativity](https://books.google.com.br/books/about/Numerical_Relativity.html?id=dxU1OEinvRUC&redir_esc=y): Section 6.2.2 (Useful to understand how to solve the Hamiltonian constraint)\n",
"\n",
"\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"1. [Step 1](#initial_data) Setting up time-symmetric initial data\n",
" 1. [Step 1.a](#id_time_symmetry) Time symmetry: $\\tilde{K}_{ij}$, $\\tilde K$, $\\tilde\\beta^{i}$, and $\\tilde B^{i}$\n",
" 1. [Step 1.b](#id_sf_ic) The scalar field initial condition: $\\tilde{\\varphi}$, $\\tilde{\\Phi}$, $\\tilde{\\Pi}$\n",
" 1. [Step 1.c](#id_metric) The physical metric: $\\tilde{\\gamma}_{ij}$\n",
" 1. [Step 1.c.i](#id_conformal_metric) The conformal metric $\\bar\\gamma_{ij}$\n",
" 1. [Step 1.c.ii](#id_hamiltonian_constraint) Solving the Hamiltonian constraint\n",
" 1. [Step 1.c.ii.1](#id_tridiagonal_matrix) The tridiagonal matrix: $A$\n",
" 1. [Step 1.c.ii.2](#id_tridiagonal_rhs) The right-hand side of the linear system: $\\vec{s}$\n",
" 1. [Step 1.c.ii.3](#id_conformal_factor) The conformal factor: $\\psi$\n",
" 1. [Step 1.d](#id_lapse_function) The lapse function: $\\tilde{\\alpha}$\n",
" 1. [Step 1.e](#id_output) Outputting the initial data to file\n",
"1. [Step 2](#id_interpolation_files) Interpolating the initial data file as needed\n",
"1. [Step 3](#id_sph_to_curvilinear) Converting Spherical initial data to Curvilinear initial data\n",
"1. [Step 4](#validation) Validation of this tutorial against the [ScalarField/ScalarField_InitialData.py](../edit/ScalarField/ScalarField_InitialData.py) module\n",
"1. [Step 5](#output_to_pdf) Output this module as $\\LaTeX$-formatted PDF file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 0: Initialize Python/NRPy+ modules \\[Back to [top](#toc)\\]\n",
"$$\\label{initialize_nrpy}$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Step 0: Load all needed Python/NRPy+ modules\n",
"import os,sys,shutil # Standard Python modules for multiplatform OS-level functions\n",
"import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
"import numpy as np # NumPy: A large collection of mathematical functions for Python\n",
"from scipy.sparse import spdiags # SciPy: Sparse, tri-diagonal matrix setup function\n",
"from scipy.sparse import csc_matrix # SciPy: Sparse matrix optimization function\n",
"from scipy.sparse.linalg import spsolve # SciPy: Solver of linear systems involving sparse matrices\n",
"import outputC as outC # NRPy+: Core C code output module\n",
"import reference_metric as rfm # NRPy+: Reference metric support\n",
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"\n",
"# Step 0.a: Create the output directory\n",
"Ccodesdir = \"ScalarFieldID_validation\"\n",
"shutil.rmtree(Ccodesdir,ignore_errors=True)\n",
"cmd.mkdir(Ccodesdir)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 1: Setting up time-symmetric initial data \\[Back to [top](#toc)\\]\n",
"$$\\label{initial_data}$$\n",
"\n",
"In this section we will set up time symmetric initial data for the gravitational collapse of a massless scalar field, in spherical coordinates. Our discussion will follow closely section III.A of [Akbarian & Choptuik (2015)](https://arxiv.org/pdf/1508.01614.pdf) (henceforth A&C). We will be using a *uniform* radial sampling. All initial data quantities will be written with tildes over them, meaning that, for example, $\\tilde{\\alpha} \\equiv \\alpha(0,r)$.\n",
"\n",
"\n",
"\n",
"## Step 1.a: Time symmetry: $\\tilde{K}_{ij}$, $\\tilde K$, $\\tilde\\beta^{i}$, and $\\tilde B^{i}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{id_time_symmetry}$$\n",
"\n",
"We are here considering a spherically symmetric problem, so that $f=f(t,r)$, for every function discussed in this tutorial. The demand for time-symmetric initial data then imples that\n",
"\n",
"\\begin{align}\n",
"\\tilde K_{ij} &= 0\\ ,\\\\\n",
"\\tilde K &= 0\\ ,\\\\\n",
"\\tilde \\beta^{i} &= 0\\ ,\\\\\n",
"\\tilde B^{i} &= 0\\ .\n",
"\\end{align}\n",
"\n",
"For the scalar field, $\\varphi$, it also demands\n",
"\n",
"$$\n",
"\\partial_{t}\\varphi(0,r) = 0\\ ,\n",
"$$\n",
"\n",
"which we discuss below.\n",
"\n",
"\n",
"\n",
"## Step 1.b: The scalar field initial condition: $\\tilde{\\varphi}$, $\\tilde{\\Phi}$, $\\tilde{\\Pi}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{id_sf_ic}$$\n",
"\n",
"We will be implementing the following options for the initial profile of the scalar field\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"\\tilde{\\varphi}_{\\rm I} &= \\varphi_{0}\\exp\\left(-\\frac{r^{2}}{\\sigma^{2}}\\right)\\ ,\\\\\n",
"\\tilde{\\varphi}_{\\rm II} &= \\varphi_{0}r^{3}\\exp\\left[-\\left(\\frac{r-r_{0}}{\\sigma}\\right)^{2}\\right]\\ ,\\\\\n",
"\\tilde{\\varphi}_{\\rm III} &= \\varphi_{0}\\left\\{1 - \\tanh\\left[\\left(\\frac{r-r_{0}}{\\sigma}\\right)^{2}\\right]\\right\\}.\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"We introduce the two auxiliary fields\n",
"\n",
"$$\n",
"\\tilde\\Phi\\equiv\\partial_{r}\\tilde\\varphi\\quad \\text{and}\\quad \\Pi\\equiv-\\frac{1}{\\alpha}\\left(\\partial_{t}\\varphi - \\beta^{i}\\partial_{i}\\varphi\\right)\\ ,\n",
"$$\n",
"\n",
"of which $\\tilde\\Phi$ will only be used as an auxiliary variable for setting the initial data, but $\\Pi$ is a dynamical variable which will be evolved in time. Because we are setting time-symmetric initial data, $\\partial_{t}\\sf = 0 = \\beta^{i}$, and thus $\\tilde\\Pi=0$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:08.486700Z",
"iopub.status.busy": "2021-06-15T10:05:08.485053Z",
"iopub.status.idle": "2021-06-15T10:05:08.820870Z",
"shell.execute_reply": "2021-06-15T10:05:08.821209Z"
}
},
"outputs": [],
"source": [
"# Step 1: Setting up time-symmetric initial data\n",
"# Step 1.a: Define basic parameters\n",
"# Step 1.a.i: Domain size\n",
"RMAX = 50\n",
"\n",
"# Step 1.a.ii: Number of gridpoints in the radial direction\n",
"NR = 30000\n",
"\n",
"# Step 1.a.iii: Initial data family. Available options are:\n",
"# Gaussian_pulse, Gaussian_pulsev2, and Tanh_pulse\n",
"ID_Family = \"Gaussian_pulsev2\"\n",
"\n",
"# Step 1.a.iv: Coordinate system. Available options are:\n",
"# Spherical and SinhSpherical\n",
"CoordSystem = \"Spherical\"\n",
"\n",
"# Step 1.a.v: SinhSpherical parameters\n",
"sinhA = RMAX\n",
"sinhW = 0.1\n",
"\n",
"# Step 1.b: Set the radial array\n",
"if CoordSystem == \"Spherical\":\n",
" r = np.linspace(0,RMAX,NR+1) # Set the r array\n",
" dr = np.zeros(NR)\n",
" for i in range(NR):\n",
" dr[i] = r[1]-r[0]\n",
" r = np.delete(r-dr[0]/2,0) # Shift the vector by -dr/2 and remove the negative entry\n",
"elif CoordSystem == \"SinhSpherical\":\n",
" if sinhA is None or sinhW is None:\n",
" print(\"Error: SinhSpherical coordinates require initialization of both sinhA and sinhW\")\n",
" sys.exit(1)\n",
" else:\n",
" x = np.linspace(0,1.0,NR+1)\n",
" dx = 1.0/(NR+1)\n",
" x = np.delete(x-dx/2,0) # Shift the vector by -dx/2 and remove the negative entry\n",
" r = sinhA * np.sinh( x/sinhW ) / np.sinh( 1.0/sinhW )\n",
" dr = sinhA * np.cosh( x/sinhW ) / np.sinh( 1.0/sinhW ) * dx\n",
"else:\n",
" print(\"Error: Unknown coordinate system\")\n",
" sys.exit(1)\n",
"\n",
"# Step 1.c: Step size squared\n",
"dr2 = dr**2\n",
"\n",
"# Step 1.d: Set SymPy variables for the initial condition\n",
"phi0,rr,rr0,sigma = sp.symbols(\"phi0 rr rr0 sigma\",real=True)\n",
"\n",
"# Step 1.e: Now set the initial profile of the scalar field\n",
"if ID_Family == \"Gaussian_pulse\":\n",
" phiID = phi0 * sp.exp( -r**2/sigma**2 )\n",
"elif ID_Family == \"Gaussian_pulsev2\":\n",
" phiID = phi0 * rr**3 * sp.exp( -(rr-rr0)**2/sigma**2 )\n",
"elif ID_Family == \"Tanh_pulse\":\n",
" phiID = phi0 * ( 1 - sp.tanh( (rr-rr0)**2/sigma**2 ) )\n",
"else:\n",
" print(\"Unkown initial data family: \",ID_Family)\n",
" print(\"Available options are: Gaussian_pulse, Gaussian_pulsev2, and Tanh_pulse\")\n",
" sys.exit(1)\n",
"\n",
"# Step 1.f: Compute Phi := \\partial_{r}phi\n",
"PhiID = sp.diff(phiID,rr)\n",
"\n",
"# Step 1.g: Generate NumPy functions for phi\n",
"# and Phi from the SymPy variables.\n",
"phi = sp.lambdify((phi0,rr,rr0,sigma),phiID)\n",
"Phi = sp.lambdify((phi0,rr,rr0,sigma),PhiID)\n",
"\n",
"# Step 1.h: populating the varphi(0,r) array\n",
"phi0 = 0.1\n",
"r0 = 0\n",
"sigma = 1\n",
"ID_sf = phi(phi0,r,r0,sigma)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 1.c: The physical metric: $\\tilde{\\gamma}_{ij}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{id_metric}$$\n",
"\n",
"\n",
"\n",
"### Step 1.c.i: The conformal metric $\\bar\\gamma_{ij}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{id_conformal_metric}$$\n",
"\n",
"To set up the physical metric initial data, $\\tilde\\gamma_{ij}$, we will start by considering the conformal transformation\n",
"\n",
"$$\n",
"\\gamma_{ij} = e^{4\\phi}\\bar\\gamma_{ij}\\ ,\n",
"$$\n",
"\n",
"where $\\bar\\gamma_{ij}$ is the conformal metric and $e^{\\phi}$ is the conformal factor. We then fix the initial value of $\\bar\\gamma_{ij}$ according to eqs. (32) and (43) of [A&C](https://arxiv.org/pdf/1508.01614.pdf)\n",
"\n",
"$$\n",
"\\bar\\gamma_{ij} = \\hat\\gamma_{ij}\\ ,\n",
"$$\n",
"\n",
"where $\\hat\\gamma_{ij}$ is the *reference metric*, which is the flat metric in spherical symmetry\n",
"\n",
"$$\n",
"\\hat\\gamma_{ij}\n",
"=\n",
"\\begin{pmatrix}\n",
"1 & 0 & 0\\\\\n",
"0 & r^{2} & 0\\\\\n",
"0 & 0 & r^{2}\\sin^{2}\\theta\n",
"\\end{pmatrix}\\ .\n",
"$$\n",
"\n",
"To determine the physical metric, we must then determine the conformal factor $e^{\\phi}$. This is done by solving the Hamiltonian constraint (cf. eq. (12) of [Baumgarte](https://arxiv.org/pdf/1807.10342.pdf))\n",
"\n",
"$$\n",
"\\hat\\gamma^{ij}\\hat D_{i}\\hat D_{j}\\psi = -2\\pi\\psi^{5}\\rho\\ ,\n",
"$$\n",
"\n",
"where $\\psi\\equiv e^{\\tilde\\phi}$. For a massless scalar field, we know that\n",
"\n",
"$$\n",
"T^{\\mu\\nu} = \\partial^{\\mu}\\varphi\\partial^{\\nu}\\varphi - \\frac{1}{2}g^{\\mu\\nu}\\left(\\partial^{\\lambda}\\varphi\\partial_{\\lambda}\\varphi\\right)\\ .\n",
"$$\n",
"\n",
"where $g^{\\mu\\nu}$ is the inverse of the ADM 4-metric given by eq. (2.119) of [Baumgarte & Shapiro's Numerical Relativity](https://books.google.com.br/books/about/Numerical_Relativity.html?id=dxU1OEinvRUC&redir_esc=y),\n",
"\n",
"$$\n",
"g^{\\mu\\nu}=\\begin{pmatrix}\n",
"-\\alpha^{-2} & \\alpha^{-2}\\beta^{i}\\\\\n",
"\\alpha^{-2}\\beta^{j} & \\gamma^{ij} - \\alpha^{-2}\\beta^{i}\\beta^{j}\n",
"\\end{pmatrix}\\ .\n",
"$$\n",
"\n",
"We know that (see Step 2 in [this tutorial module](Tutorial-ADM_Setting_up_massless_scalarfield_Tmunu.ipynb) for the details)\n",
"\n",
"\\begin{align}\n",
"\\partial^{t}\\varphi &= \\alpha^{-1}\\Pi\\ ,\\\\\n",
"\\partial^{\\lambda}\\varphi\\partial_{\\lambda}\\varphi &= -\\Pi^{2} + \\gamma^{ij}\\partial_{i}\\varphi\\partial_{j}\\varphi\\ .\n",
"\\end{align}\n",
"\n",
"The tt-component of the energy-momentum tensor at the initial time is then given by (we will ommit the \"tildes\" below to avoid cluttering the equation, but keep in mind that all quantities are considered at $t=0$)\n",
"\n",
"\\begin{align}\n",
" T^{tt} &= \\left(\\partial^{t}\\varphi\\right)^{2} - \\frac{1}{2} g^{tt}\\left(\\partial^{\\lambda}\\varphi\\partial_{\\lambda}\\varphi\\right)\\nonumber\\\\\n",
"&= \\left(\\frac{\\Pi}{\\alpha}\\right)^{2} - \\frac{1}{2}\\left(-\\frac{1}{\\alpha^{2}}\\right)\\left(-\\Pi^{2} + \\gamma^{ij}\\partial_{i}\\varphi\\partial_{j}\\varphi\\right)\\nonumber\\\\\n",
"&= \\frac{\\gamma^{ij}\\partial_{i}\\varphi\\partial_{j}\\varphi}{2\\alpha^{2}}\\nonumber\\\\\n",
"&= \\frac{e^{-4\\phi}\\bar\\gamma^{ij}\\partial_{i}\\varphi\\partial_{j}\\varphi}{2\\alpha^{2}}\\nonumber\\\\\n",
"&= \\frac{e^{-4\\phi}\\hat\\gamma^{ij}\\partial_{i}\\varphi\\partial_{j}\\varphi}{2\\alpha^{2}}\\nonumber\\\\\n",
"&= \\frac{e^{-4\\phi}\\hat\\gamma^{rr}\\partial_{r}\\varphi\\partial_{r}\\varphi}{2\\alpha^{2}}\\nonumber\\\\\n",
"&= \\frac{e^{-4\\phi}\\Phi^{2}}{2\\alpha^{2}}\\nonumber\\\\\n",
"\\end{align}\n",
"\n",
"By remembering the definition of the normal vector $n_{\\mu} = (-\\alpha,0,0,0)$ (eq. (2.117) of [B&S](https://books.google.com.br/books/about/Numerical_Relativity.html?id=dxU1OEinvRUC&redir_esc=y)), we can then evaluate the energy density $\\rho$ given by eq. (24) of [A&C](https://arxiv.org/pdf/1508.01614.pdf)\n",
"\n",
"$$\n",
"\\tilde\\rho = \\tilde n_{\\mu}\\tilde n_{\\nu}\\tilde T^{\\mu\\nu} = \\frac{e^{-4\\tilde\\phi}}{2}\\tilde\\Phi^{2}\\ .\n",
"$$\n",
"\n",
"Plugging this result in the Hamiltonian constraint, remembering that $\\psi\\equiv e^{\\tilde\\phi}$, we have\n",
"\n",
"$$\n",
"\\partial^{2}_{r}\\psi + \\frac{2}{r}\\partial_{r}\\psi + \\pi\\psi\\Phi^{2} = 0\\ .\n",
"$$\n",
"\n",
"This is a linear elliptic equation which will solve using the procedure described in detail in section 6.2.2 of [B&S](https://books.google.com.br/books/about/Numerical_Relativity.html?id=dxU1OEinvRUC&redir_esc=y).\n",
"\n",
"\n",
"\n",
"### Step 1.c.ii: Solving the Hamiltonian constraint \\[Back to [top](#toc)\\]\n",
"$$\\label{id_hamiltonian_constraint}$$\n",
"\n",
"We will discretize the Hamiltonian constraint using [second-order accurate finite differences](https://en.wikipedia.org/wiki/Finite_difference_coefficient). We get\n",
"\n",
"$$\n",
"\\frac{\\psi_{i+1} - 2\\psi_{i} + \\psi_{i-1}}{\\Delta r^{2}} + \\frac{2}{r_{i}}\\left(\\frac{\\psi_{i+1}-\\psi_{i-1}}{2\\Delta r}\\right) + \\pi\\psi_{i}\\Phi^{2}_{i} = 0\\ ,\n",
"$$\n",
"\n",
"or, by multiplying the entire equation by $\\Delta r^{2}$ and then grouping the coefficients of each $\\psi_{j}$:\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 choose to set up a grid that is cell-centered, with:\n",
"\n",
"$$\n",
"r_{i} = \\left(i-\\frac{1}{2}\\right)\\Delta r\\ ,\n",
"$$\n",
"\n",
"so that $r_{0} = - \\frac{\\Delta r}{2}$. This is a two-point boundary value problem, which we solve using the same strategy as [A&C](https://arxiv.org/pdf/1508.01614.pdf), described in eqs. (48)-(50):\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",
"In terms of our grid structure, the first boundary condition (regularity at the origin) is written to second-order in $\\Delta r$ as:\n",
"\n",
"$$\n",
"\\left.\\partial_{r}\\psi\\right|_{r=0} = \\frac{\\psi_{1} - \\psi_{0}}{\\Delta r} = 0 \\Rightarrow \\psi_{0} = \\psi_{1}\\ .\n",
"$$\n",
"\n",
"The second boundary condition (asymptotic flatness) can be interpreted as\n",
"\n",
"$$\n",
"\\psi_{N} = 1 + \\frac{C}{r_{N}}\\ (r_{N}\\gg1)\\ ,\n",
"$$\n",
"\n",
"which then implies\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 then 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",
"Substituting the boundary conditions at the boxed equations above, we end up with\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\n",
"\n",
"#### Step 1.c.ii.1: The tridiagonal matrix: $A$ \\[Back to [top](#toc)\\]\n",
"$$\\label{id_tridiagonal_matrix}$$\n",
"\n",
"We now start solving the tridiagonal linear system. We start by implementing the tridiagonal matrix $A$ defined above. We break down it down by implementing each diagonal into an array. We start by looking at the main diagonal:\n",
"\n",
"$$\n",
"{\\rm diag}_{\\rm main}\n",
"=\n",
"\\begin{pmatrix}\n",
"\\left(\\pi\\Delta r^{2}\\Phi^{2}_{1} - 1 - \\frac{\\Delta r}{r_{1}}\\right)\\\\\n",
"\\left(\\pi\\Delta r^{2}\\Phi_{2}^{2}-2\\right)\\\\\n",
"\\vdots\\\\\n",
"\\left(\\pi\\Delta r^{2}\\Phi_{i}^{2}-2\\right)\\\\\n",
"\\vdots\\\\\n",
"\\left(\\pi\\Delta r^{2}\\Phi_{N-1}^{2}-2\\right)\\\\\n",
"\\left[\\pi\\Delta r^{2}\\Phi^{2}_{N} - 2 - \\frac{2\\Delta r}{r_{N}}\\left(1+\\frac{\\Delta r}{r_{N}}\\right)\\right]\\\\\n",
"\\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
"\\left(\\pi\\Delta r^{2}\\Phi^{2}_{1} - 2\\right)\\\\\n",
"\\left(\\pi\\Delta r^{2}\\Phi_{2}^{2} - 2\\right)\\\\\n",
"\\vdots\\\\\n",
"\\left(\\pi\\Delta r^{2}\\Phi_{i}^{2} - 2\\right)\\\\\n",
"\\vdots\\\\\n",
"\\left(\\pi\\Delta r^{2}\\Phi_{N-1}^{2}-2\\right)\\\\\n",
"\\left(\\pi\\Delta r^{2}\\Phi^{2}_{N} - 2\\right)\\\\\n",
"\\end{pmatrix}\n",
"+\n",
"\\left.\\begin{pmatrix}\n",
"1 - \\frac{\\Delta r}{r_{1}}\\\\\n",
"0\\\\\n",
"\\vdots\\\\\n",
"0\\\\\n",
"\\vdots\\\\\n",
"0\\\\\n",
"- \\frac{2\\Delta r}{r_{N}}\\left(1+\\frac{\\Delta r}{r_{N}}\\right)\n",
"\\end{pmatrix}\\quad \\right\\}\\text{N elements}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:08.823299Z",
"iopub.status.busy": "2021-06-15T10:05:08.822859Z",
"iopub.status.idle": "2021-06-15T10:05:08.828549Z",
"shell.execute_reply": "2021-06-15T10:05:08.828895Z"
}
},
"outputs": [],
"source": [
"# Set the main diagonal\n",
"main_diag = np.pi * dr2 * Phi(phi0,r,r0,sigma)**2 - 2\n",
"\n",
"# Update the first element of the main diagonal\n",
"main_diag[0] += 1 - dr[0]/r[0]\n",
"\n",
"# Update the last element of the main diagonal\n",
"main_diag[NR-1] += - (2 * dr[NR-1] / r[NR-1])*(1 + dr[NR-1] / r[NR-1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we look at the upper diagonal of the A matrix:\n",
"\n",
"$$\n",
"{\\rm diag}_{\\rm upper}\n",
"=\n",
"\\left.\\begin{pmatrix}\n",
"1+\\frac{\\Delta r}{r_{1}}\\\\\n",
"1+\\frac{\\Delta r}{r_{2}}\\\\\n",
"\\vdots\\\\\n",
"1+\\frac{\\Delta r}{r_{i}}\\\\\n",
"\\vdots\\\\\n",
"1+\\frac{\\Delta r}{r_{N-2}}\\\\\n",
"1+\\frac{\\Delta r}{r_{N-1}}\n",
"\\end{pmatrix}\\quad\\right\\}\\text{N-1 elements}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:08.830461Z",
"iopub.status.busy": "2021-06-15T10:05:08.830102Z",
"iopub.status.idle": "2021-06-15T10:05:08.832191Z",
"shell.execute_reply": "2021-06-15T10:05:08.832537Z"
}
},
"outputs": [],
"source": [
"# Set the upper diagonal, ignoring the last point in the r array\n",
"upper_diag = np.zeros(NR)\n",
"upper_diag[1:] = 1 + dr[:-1]/r[:-1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we look at the lower diagonal of the A matrix:\n",
"\n",
"$$\n",
"{\\rm diag}_{\\rm lower}\n",
"=\n",
"\\left.\\begin{pmatrix}\n",
"1-\\frac{\\Delta r}{r_{2}}\\\\\n",
"1-\\frac{\\Delta r}{r_{3}}\\\\\n",
"\\vdots\\\\\n",
"1-\\frac{\\Delta r}{r_{i+1}}\\\\\n",
"\\vdots\\\\\n",
"1-\\frac{\\Delta r}{r_{N-1}}\\\\\n",
"2\n",
"\\end{pmatrix}\\quad\\right\\}\\text{N-1 elements}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:08.833995Z",
"iopub.status.busy": "2021-06-15T10:05:08.833640Z",
"iopub.status.idle": "2021-06-15T10:05:08.835879Z",
"shell.execute_reply": "2021-06-15T10:05:08.836212Z"
}
},
"outputs": [],
"source": [
"# Set the lower diagonal, start counting the r array at the second element\n",
"lower_diag = np.zeros(NR)\n",
"lower_diag[:-1] = 1 - dr[1:]/r[1:]\n",
"\n",
"# Change the last term in the lower diagonal to its correct value\n",
"lower_diag[NR-2] = 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we construct the tridiagonal matrix by adding the three diagonals, while shifting the upper and lower diagonals to the right and left, respectively. Because A is a sparse matrix, we will also use scipy to solve the linear system faster."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:08.837659Z",
"iopub.status.busy": "2021-06-15T10:05:08.837306Z",
"iopub.status.idle": "2021-06-15T10:05:09.418611Z",
"shell.execute_reply": "2021-06-15T10:05:09.418110Z"
}
},
"outputs": [],
"source": [
"!pip install scipy >/dev/null"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:09.422813Z",
"iopub.status.busy": "2021-06-15T10:05:09.422294Z",
"iopub.status.idle": "2021-06-15T10:05:09.427756Z",
"shell.execute_reply": "2021-06-15T10:05:09.427226Z"
}
},
"outputs": [],
"source": [
"# Set the sparse matrix A by adding up the three diagonals\n",
"A = spdiags([main_diag,upper_diag,lower_diag],[0,1,-1],NR,NR)\n",
"\n",
"# Then compress the sparse matrix A column wise, so that SciPy can invert it later\n",
"A = csc_matrix(A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"#### Step 1.c.ii.2 The right-hand side of the linear system: $\\vec{s}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{id_tridiagonal_rhs}$$\n",
"\n",
"We now focus our attention to the implementation of the $\\vec{s}$ vector:\n",
"\n",
"$$\n",
"\\vec{s} = \n",
"\\begin{pmatrix}\n",
"0\\\\\n",
"0\\\\\n",
"\\vdots\\\\\n",
"0\\\\\n",
"\\vdots\\\\\n",
"0\\\\\n",
"-\\frac{2\\Delta r}{r_{N}}\\left(1+\\frac{\\Delta r}{r_{N}}\\right)\n",
"\\end{pmatrix}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:09.430481Z",
"iopub.status.busy": "2021-06-15T10:05:09.429977Z",
"iopub.status.idle": "2021-06-15T10:05:09.432341Z",
"shell.execute_reply": "2021-06-15T10:05:09.431827Z"
}
},
"outputs": [],
"source": [
"# Set up the right-hand side of the linear system: s\n",
"s = np.zeros(NR)\n",
"\n",
"# Update the last entry of the vector s\n",
"s[NR-1] = - (2 * dr[NR-1] / r[NR-1])*(1 + dr[NR-1] / r[NR-1])\n",
"\n",
"# Compress the vector s column-wise\n",
"s = csc_matrix(s)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"#### Step 1.c.ii.3 The conformal factor: $\\psi$ \\[Back to [top](#toc)\\]\n",
"$$\\label{id_conformal_factor}$$\n",
"\n",
"We now use scipy to solve the sparse linear system of equations and determine the conformal factor $\\psi$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:09.434659Z",
"iopub.status.busy": "2021-06-15T10:05:09.434155Z",
"iopub.status.idle": "2021-06-15T10:05:09.454515Z",
"shell.execute_reply": "2021-06-15T10:05:09.454110Z"
}
},
"outputs": [],
"source": [
"# Solve the sparse linear system using scipy\n",
"# https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.spsolve.html\n",
"psi = spsolve(A, s.T)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We then show useful plots of the conformal factor $\\psi$ and of the *evolved conformal factors*\n",
"\n",
"\\begin{align}\n",
"\\phi &= \\log\\psi\\ ,\\\\\n",
"W &= \\psi^{-2}\\ ,\\\\\n",
"\\chi &= \\psi^{-4}\\ .\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:09.460946Z",
"iopub.status.busy": "2021-06-15T10:05:09.460411Z",
"iopub.status.idle": "2021-06-15T10:05:10.017251Z",
"shell.execute_reply": "2021-06-15T10:05:10.017683Z"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABLAAAAMgCAYAAAAz4JsCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAA9hAAAPYQGoP6dpAADL30lEQVR4nOzdd7gdVbn48e+b5CQkIQmhJlRDR+lVsCEIIiqiggV/ImBD4SpXEHvDe0XxiiDYQAQsIGDBRgkIKAhSBektoSdAKElIPUne3x8zO9k5nLZz9sneZ5/v53nm2XvPrJl5ZwWyVt5ZsyYyE0mSJEmSJKlZDWl0AJIkSZIkSVJ3TGBJkiRJkiSpqZnAkiRJkiRJUlMzgSVJkiRJkqSmZgJLkiRJkiRJTc0EliRJkiRJkpqaCSxJkiRJkiQ1NRNYkiRJkiRJamomsCRJkiRJktTUTGBJkiRJkiSpqZnAkiRJkiRJUlMzgSVJkiRJkqSmZgJLkiRJkiRJTc0EliRJkiRJkpqaCSxJkiRJkiQ1NRNYkiRJkiRJamomsCRJkiRJktTUTGBJkiRJkiSpqZnAkiRJkiRJUlMzgSVJkiRJkqSmZgJLkiRJkiRJTc0EliRJkiRJkpqaCSxJkiRJkiQ1NRNYkiRJkiRJamomsCRJkiRJktTUTGBJkiRJkiSpqZnAkiRJkiRJUlMzgSVJkiRJkqSmZgJLkiRJkiRJTc0EliRJkiRJkpqaCSxJkiRJkiQ1NRNYkiRJkiRJamomsCRJkiRJktTUTGBJkiRJkiSpqZnAkiRJkiRJUlMzgSVJkiRJkqSmZgJLkiRJkiRJTc0EliRJkiRJkpqaCSxJkiRJkiQ1NRNYkiRJkiRJamomsCRJkiRJktTUTGBJkiRJkiSpqZnAkiRJkiRJUlMzgSVJkiRJkqSmZgJLkiRJkiRJTc0EliRJkiRJkpqaCSxJkiRJkiQ1NRNYkiRJkiRJamomsCRJkiRJktTUTGBJkiRJkiSpqZnAkiRJkiRJUlMzgSVJkiRJkqSmZgJLkiRJkiRJTc0EliRJkiRJkpqaCSxJkiRJkiQ1NRNYkiRJkiRJamomsCRJkiRJktTUTGBJkiRJkiSpqZnAkiRJkiRJUlMzgSVJkiRJkqSmZgJLkiRJkiRJTc0EliRJkiRJkpqaCSxJkiRJkiQ1NRNYkiRJkiRJamomsCRJkiRJktTUTGBJkiRJkiSpqZnAkiRJkiRJUlMzgSVJkiRJkqSmZgJLkiRJkiRJTc0EliRJkiRJkpqaCSxJkiRJkiQ1NRNYkiRJkiRJamomsCRJkiRJktTUTGBJanoRsUtEXB8RcyIiI2L7RsdUi4g4rIz7Fb0oO6CvtV4i4viIuC8iVmo7FRFHRsRjETFiZZ5XklamWtqlfjj31yMiV/Z5e6sV2mH7HS8XEbdGxB+62NaQPkd5bvsdUg1MYEnqlYjYJCJ+GhFTImJ+RMyKiH9GxKcjYmQ/nrcNuAhYHfhv4IPAo/11vkZaGdcaEXuU/3hYrZ7HraeIGAt8DvhOZi6pWj8iIr4TEU9FxLyIuDEi9qnz6c8BhgMfr/NxJellqhINXS2vbnSMg8lg6nPA4Ol3REQAWwL3dLKt0z5Huc1+h9RkhjU6AEnNLyLeStHBWQD8AriLorF9LfBd4FXAx/rp9JsAGwEfzcyf9dM5msXKuNY9gK9RdJhe7Kdz9NURFO3T+R3WnwMcBJwCPAgcBlwSEW/MzOvqceLMnB8R5wKfiYjTMrNpRwlIailfBaZ2sv6hlR3IIDeY+hwwePodrwBG0UkCi677HGC/Q2o6JrAkdSsiJgG/obgjt1dmTqva/MOI2BR4az+GsHb5+WK9DhgRozNzTr2OV0d1v9aVpc51ejjwp8ycX3X8XYH3AZ/NzP8r11WSqSdRdJD7LCJWAS4EjgfeCFxVj+NKUg8uzcxbGh2EBlWfAwZov2MF6vSV5WdnCayX9TnKc/R7vyMiVinPa79D6iUfIZTUk+OBVYEPd0heAZCZD2XmqZXfEbFDRFxaPmL4UkT8reMjEJX5LyJi04g4JyJejIiZEXF2RIyqKncO8Pfy50XlPtes4LleGRHnRcQLwHUdtm0eEb8qY3g2Ir4ZhQ0i4o/l8adHxLEdjr1RRPwoIu4vh5Y/FxEXxQrMKdLdtdZynohYLyLOKoe7L4iIqRHx44gYHhFfpxgxBzC16hGVV1Tt36c67eLaojzW98rfm5f7H1z+PrT8vVmZMN0WuLLDYQ4CFgNnVFaUnb6zgN0jYoNuK7jzuK6IYt6P10XE3yNiHnBqZt4KPA+8o9ZjSlK9RcRB5d+Rb+hk28fLbVtXrevx7/G+nKNsZ34eEU+X7czdEXFEF8d9bUTcHMXUAw9HRE2PSXXXptVyvdGLfkd37fAKnKep+xzdXW+t5+nuzyh66Hf0tU57cY3vjIhbgcrcV9dGxK8jYly5vas+B9S53xFd9DnK49rvkHrJEViSevJ2YEpmXt9TwYh4FXAtMIvi7lQ7xTP910TEGzLzxg67XEjxyMQXgB2BjwDPUMxFAPBT4Engi8APgJuBp1fwXBdRDP/+IhAdtl0A3At8nmI02ZcpOhIfp7gT9jngA8D/RcTNmfmPcr9dKO7A/QZ4gmKI+ifKGF6ZmXN7qrMqXV5rb88TEesCNwGrUXS47gPWo+iEjQJ+D2wOvJ9irosZ5fGfLfevZ51W2xgYDdxZ/t6m/Kz+PQ94uIwN4LYOx9gBeCAzZ3VYf1P5uT3weDcxdGZbirvOFwNnAuex7HGd24DX1Hg8SVpR4yJizQ7rMjOfA/4KvAS8h2UJh4r3Andn5l2wQn+PV/TqHBGxDvAvIIHTKdqPtwBnRcTYzDylslNEbANMLst8neLfHd9gWdvWrV60aQvr3O8YTH0OurneXp+nL/2OfuxzVGL7bHnc84E2YCxwWXmOAA5h2Siqjn0OqH+/o7s+RyUG+x1STzLTxcXFpdOForFP4OJelv8DxTxZG1etm0jROfl71bqvl8c9q8P+vwdmdFi3Z1n2oD6e67xO4q1s+2nVuqEUHZIlwOeq1q8GzAXOqVo3spNjvro85ger1h1WrntFD/XX1bX29jznUtwt3LmT8lF+HtdVLPWo0y6u6x1l+Z2q9p8PDC1/XwrcXH7/Zll21Q7HuAv4WyfHfmVZ/uM1/re9drnfbGDLTrb/FJjbn/9/ubi4uFS1D50t86vKnUeRXBhatW5C+Xf+V6rW9fbv8Ze1S705B/Az4ClgjQ7XcT7FP85HdohlHrBh1bqtgEUUybme6qY3bVpd+x20QJ+jqz/fLur4Zddb43lWuN9Rjzrt5rp2Kev0u+Xv+yv7UyRV2ykSbJ32Ocpydet30EOfoyxjv8PFpReLjxBK6s7Y8nN2TwUjYiiwL0Wya0plfRaPHZ4HvDaKN71U+0mH39cCa3RSrj/OVW3pxKWZuRi4heLu3FlV61+k6ABtXLVuXlVMbRGxBsXdtBcp7uzWRW/OE8Wrnw8E/pydzKOSmdndOfqhTqttTdGRvKfq971lXVd+V0ZjrQEsysyXOhxjJEVHt6P5VdtrsW35+a3MvK+T7S8AI6PqkVZJ6kdHAft0WN5Stf0Cin8E71m17iCK6UAugBX+e7xat+eIiADeDfy5OF2sWVmAy4FxLGuThgJvLmN5rCqWe8uy3epNm7ay+h32OTo/T1/6Hf3c54BiFNuzwFejeFP2psAd5bZ/UowGXJuu+xxQ335HT30OsN8h9YoJLEndqQybHtOLsmtR3M26v5Nt91L8fdNxvoDHOvx+ofwc3w/nmtrN8TrGMZPizveMTtYvjS0iRkbECRHxOEUnZwZFh2k1io58XfTyPGtRJBzvWsHT1LtOq20DPFzVKd6GMmFVzkOxPssSWF2ZB4zoZP0qVdtrUXmM8YIutlceT+g28SdJdXJTZl7ZYbm6avtlFG3Qe6vWvRe4PTMfKH+vyN/j1Xo6x1oU7c7HKNqg6uXssnxlUvC1KP6B/2An5+ksvo5606atrH6HfY7Oz9OXfke/9TkiYhiwH8WLEeZR3CQbAvynLDK6/Hyhk92r1bPf0VOfA+x3SL3iHFiSupSZsyLiKYrGvz8s7mJ9j3MbrIDuOhqdxdGb2E6jeHvNKcANFJ3NpJg3op43CFbWeWrV287bqyg7uFG85W8Tlt1prvy3VUlgPQcMi4gxmVk98m8axbwaHU0sP5/qbdClbYFp1Xd+OxhPMZS/1sSYJNVdZi6IiIuBd0bEJ4F1KObL+eJKPEelvfkVxaNjnflPF+ubxcrqd9jnqL/etsebsvy8m5XRT5URWNsDj2bmzIjoqs8B9e139NTnAPsdUq+YwJLUk78AH4uI3TPzhm7KPUsxX8MWnWzbkuIRslon2W6Gc3XnIODczFz6pqAyQbNaA87zLMWIuZ6SjV3d2evPOt2Y4jENKOaOGMqyjuWbymNXHj+oDK2fxPL/ELodeGM5SXD1hKq7VW2vxbYs68x2ZhLFXWBJahYXAB8C9qaYSypYfkRHPf4e7+4cz1JMKTA0Mzt7a1u1ZykSDpt1sq2z+Drbv6c2bWX1BexzdH6evvQ7+rNOK6PW5pSf21HMc/ZU+bjrG1j2OGJXfQ6ob7+jpz5HJQb7HVIPfIRQUk9OougE/Kx8+9ByImKTiPh0OY/DZOAdUfWq5XKfQ4Dr8uVvclkhK/NcPVjMy+/a/hdFgmalniczl1C82ebtEbFzxwOUc5fAsg7datXb+7lOhwGVt2tVOrp3lY8PHg5cWc73AcXdXoCO1/Bbiuv9WFVsI8r9b8zMXnd0y7k3Xkn3nckdgR7fvClJK9GVFG+re2+53JSZSx+rqtPf412eozz+74B3R8TLkhYRsVaHWC4HDoyIDavKbEUxN1a3etOmray+gH2Ozs/Tl35HP9fpk+Xn7uVndfLo+xTJsVPK3131OaBO/Y5e9jnAfofUK47AktStzHw4Ig6hfO1zRPyC4nGw4RSvHz4YOKcs/mWKiWevi4gfUbxp6OMUcwgcX+fQVua5uvIX4IMRMZNigvLdKUYUPdeg83yRYlLUv0fEGRR38iZS/Bm9lmIC1lvLsv8bEb+heBPPnzNzDv1Xp/+k+EfMpyjmgZhLkcj6OsV8KQdWCmbmlIi4q7y+n1etvzEiLgJOjIi1KSaU/RDF670/3PGEEZEUbzHas5N4NqOYw6LTzmRE7ASsDvyxtsuUpBX2lojYspP111ceO8rM9oj4PfA+ikekjuukfJ/+Hu/FOT4PvBG4MSLOpGiTVqf4x/ebyu8VX6OYi+jaMpZhFImQu1n2WFd3etOmray+gH2Ozs+zwv0O+qlOM/OxiLimvIanKUZg3RURf6Z4McKHqpKynfY5ym316nd02+co97XfIfWSCSxJPcrMP0XEtsBngXcAn6CY2PM/wLHAmWW5uyPidcCJwBcoRnneCPy/zLyxzjGttHN149MUdyo/QNE5+SdFJ6jHNyz1x3ky88mI2I3itdAfoJhc9UngUoqkEZl5c0R8BTiS4h8WQyiGrc/pxzo9kuKu/alV6y4pr+O1mXl7h/I/B06IiJEd5oI4tLy2D1I8IvAf4G2Z+Y/qnSNi1fLrtC7iqUym2tVcLQdTTLJ7VVcXJEl1dkIX6w8HqufNuQD4CMVjWRd2LFynv8e7PEdmPh0RuwJfBd4FfJIisXE3xZvfqsv+JyLeDJxcXt8TFEmtifQigdXLNm2l9AXsc3R+nr70O/q5Tt9P0Tf9r/IadgVuBvbNzI5te1d9DqhPv6OnPgfY75B6LXp4s7okSX1WPkqwK3AF8A/gk9WvVu9QdhzFP9iOz8yzOivTw7n2p7iDvF1m9vR2w477jgAeAb6dmaf2UFySJDWpiHgbxWivHTq5WVYp06c+R3kM+x3SSuIcWJKkfpfF3ZL7gTHAZV0lr8qyMynmXvtsRKxIO/VG4De1diJLh1M83vCTngpKkqSmtiXFSML7uypQhz4H2O+QVhpHYEmSVoqIeC3F2wj3zMy/NzoeSZLUuiLiZ8CbMvMVjY5FUn04AkuStLJU5oFYkTuUkiRJtdgSuK/RQUiqH0dgSZIkSZIkqak5AkuSJEmSJElNzQSWJEmSJEmSmtqwRgeg1hMRAawLzG50LJIkVRkDPJXOn9BntvWSpCZlW9/CTGCpP6wLPNHoICRJ6sT6wJONDqIF2NZLkpqVbX2LMoGl/jAb4PHHH2fs2LF9OlB7ezuTJ09m3333pa2trS7BtTrrrHbWWe2ss9pZZ7WrZ53NmjWLDTbYABwxVC+29Q1kndXOOquddVY766x2tvWqhQks9ZuxY8fWpVM7atQoxo4dayPQS9ZZ7ayz2llntbPOamedNT/b+sawzmpnndXOOquddVY760y1cBJ3SZIkSZIkNTUTWJIkSQ0UEUdFxCMRMT8iboyIXXsof3BE3FeWvzMi9u+w/V0RMTkinouIjIjtO2xfPSJOi4j7I2JeRDwWET+IiHH9cHmSJEl1YQJLkiSpQSLivcDJwDeAHYE7gMsjYu0uyu8BnA+cBewAXAxcHBFbVxUbDVwHfK6L065bLscBWwOHAfuVx5QkSWpKJrCaVES8PiL+HBFPlXdPD+zFPntGxG0RsSAiHoqIw2o5ZkS0RcR3yru5c8pyv4iIdet6cZIkqeIzwJmZeXZm3gMcCcwFjuii/KeByzLzu5l5b2Z+BbgNOLpSIDN/mZknAFd2doDMvCsz352Zf87MhzPzKuBLwNsjwvlRJUlSUzKB1bxGU9yFPao3hSNiEvBX4Gpge+AU4GcR8eYajjmK4u7vN8vPdwFbAH+qOXpJktStiBgO7ERVoikzl5S/d+9it915eWLq8m7K99Y4YFZmLurjcSRJkvqFd9maVGZeClwKEBG92eVIYGpmHlv+vjciXgv8N0XHtsdjZuZMYJ/qdRFxNHBTRGyYmY+t0MVIkqTOrAkMBZ7usP5pYMsu9pnQRfkJKxpERKwJfAU4o5syI4ARVavGQPH2qPb29hU9NZVjVH+qZ9ZZ7ayz2llntbPOalfPOrPeW58JrNbR1R3ZU/p43HFAAi92VcBObXOxzmpnndXOOquddVY7O7X9LyLGUozgvgf4ejdFvwB8rePKyZMnM2rUqLrEcsUVV9TlOIOJdVY766x21lntrLPa1aPO5s6dW4dI1MxMYLWOru7Ijo2IkZk5r9YDRsQqwHeA8zNzVjdF7dQ2IeusdtZZ7ayz2llntWvhTu0MYDGwTof16wDTu9hneo3luxQRY4DLgNnAOzOzuyzfiRSTzVeMAZ7Yd999GTt2bK2nXk57eztXXHEF++yzD21tbX061mBhndXOOquddVY766x29ayzWbO6+yerWoEJLHUqItqAC4EAPtFD8X7r1H7rW9/itttu41WvehXjx4/nLW95C1tssUWfjtnqbDhrZ53VzjqrnXVWu1bv1Gbmwoi4Fdib4m2CRMSQ8vfpXex2Q7n9lKp1+5Tre60ceXU5sAA4IDPn9xDrgrJsZX8A2tra+vxnE3/8I1ufey4jbriBoePGwa67wl57Qe+mUBjU6lH/g411VjvrrHbWWe3qUWfWeeszgdU6urojO6vW0VdVyauNgL16GH3Vr53aCy64gHvvvZc//amYR/6LX/wiF1xwAe9617v6dNzBwIazdtZZ7ayz2llntWvxTu3JwLkRcQtwE3AMxUtXzgaIiF8AT2bmF8rypwJ/j4hjKR79ex+wM/CxygEjYnVgQ6DyFuEtyrZ5emZOL5NXkyle3vL/KEZrV+44PZuZi/vpWjsV117LJn/5C/zlL8tWHn00nHbaygxDkiQ1ORNYreMGYP8O61bkjmwlebUZ8MbMfK4+4a2Yj370o1x77bVMmDCBu+66i2uvvZaPfOQjvPGNb2T8+PGNDE2SpD7LzAsiYi3gBIrpAG4H9svMyrQAGwJLqspfHxGHAP8DfAt4EDgwM++qOuwBlAmw0m/Kz29QzHO1I7Bbue6hDiFNAh7p00XVKPfemweefJJNJk5k6DPPwEUXwemnwwEHwD779HwASZI0KJjAalIRsSqwadWqSRGxPfB8Zj4WEScC62XmoeX2nwBHR8RJwM+BvYD3AG+t4ZhtwG8pOrZvA4ZGROWtRs9n5sJ6X2dPjj76aDbeeGP2339/IoLtttuOe+65h1/+8pd86lOfWtnhSJJUd5l5Ol08MpiZe3ay7iLgom6Odw5wTjfbr6GYIqAp5Fvewr2ZTNp/f4a2tcE66xSjr045xQSWJElaakijA1CXdgb+XS5QPGLwb4o7tAATKe7KApCZUymSVfsAdwDHAh/JzMtrOOZ6FHdt16e4Azytatmjble2goYNG8aRRx4JwIUXXtjgaCRJUr84+uji8/LLYebMxsYiSZKahgmsJpWZ12RmdLIcVm4/rONd2XKfHTJzRGZuUt6BreWYj3SxPcq7tQ33tre9DYB//etfTTkhryRJ6qPNNy+WxYvhmmsaHY0kSWoSJrA0oEyaNImNN96YxYsXc/311zc6HEmS1B/23rv4NIElSZJKJrA04Oy2WzHv7L///e8eSkqSpAGpbOu57bbGxiFJkpqGCSwNODvssANgAkuSpJZVtvXcfjssWdJtUUmSNDiYwNKAs+OOOwJwm3dlJUlqTVttBSNGwKxZMGVKo6ORJElNwASWBpxtttkGgClTprBgwYIGRyNJkuqurQ223LL4ft99jY1FkiQ1BRNYGnDWWmstxo0bR2by8MMPNzocSZLUH7bYovi8//7GxiFJkpqCCSwNOBHB5ptvDsADDzzQ4GgkSVK/KNt6bOslSRImsDRAVRJY93tXVpKk1mQCS5IkVTGBpQHJEViSJLU4HyGUJElVTGBpQNpss80AE1iSJLWsygisadNg9uzGxiJJkhrOBJYGpEoC66GHHmpwJJIkqV+sthqsuWbx3Ze2SJI06JnA0oC0ySabADB9+nTmzJnT4GgkSVK/2HTT4tMbVpIkDXomsDQgjR8/njXWWAOAh70rK0lSazKBJUmSSiawNGBVRmH5GKEkSS3KBJYkSSqZwNKAtWnZqTWBJUlSizKBJUmSSiawNGCZwJIkqcWZwJIkSSUTWBqwTGBJktTiKgmsJ5+EuXMbG4skSWooE1gasExgSZLU4lZfHVZbrfg+ZUpDQ5EkSY1lAksDViWB9cQTTzB//vwGRyNJkuouwscIJUkSYAJLA9iaa67J2LFjyUymTp3a6HAkSVJ/MIElSZIwgaUBLCLYZJNNAB8jlCSpZZnAkiRJmMDSAOc8WJIktTgTWJIkCRNYGuBMYEmS1OJMYEmSJExgaYAzgSVJUourJLAeewwWLGhsLJIkqWFMYGlAqySwHn744QZHIkmS+sXaa8Oqq0Im+NIWSZIGLRNYGtAqCaxHHnmE9vb2BkcjSZLqLsLHCCVJkgksDWwTJ05k5MiRLF68mEcffbTR4UiSpP5gAkuSpEHPBJYGtIhgk002AZwHS5KklmUCS5KkQc8ElgY8J3KXJKnFmcCSJGnQM4GlAc8EliRJLc4EliRJg54JLA14volQkqQWV0lgPfII+NIWSZIGJRNYGvAcgSVJGsgi4qiIeCQi5kfEjRGxaw/lD46I+8ryd0bE/h22vysiJkfEcxGREbF9J8dYJSJ+WJZ5KSJ+FxHr1PnS6mfiRBg5EhYvBl/aIknSoGQCq0lFxOsj4s8R8VTZ+TywF/vsGRG3RcSCiHgoIg6r9ZhROCEipkXEvIi4MiI2q9uF9YPqEVjt3pWVJA0gEfFe4GTgG8COwB3A5RGxdhfl9wDOB84CdgAuBi6OiK2rio0GrgM+182pvw+8HTgYeAOwLvD7vlxLvxoyZNkorHvvbWwskiSpIUxgNa/RFJ3Yo3pTOCImAX8Frga2B04BfhYRb67xmMcDnwKOBHYD5lB0pFepLfyVZ8MNN2TcuHG0t7dz9913NzocSZJq8RngzMw8OzPvoWh/5wJHdFH+08BlmfndzLw3M78C3AYcXSmQmb/MzBOAKzs7QESMAz4MfCYzr8rMW4HDgT0i4tV1u7J622GH4vO22xobhyRJaohhjQ5AncvMS4FLASKiN7scCUzNzGPL3/dGxGuB/wYu780xo1h5DPA/mfnHct2hwNPAgcBvVvR6+lNEsOOOO3L11Vdz6623sv322zc6JEmSehQRw4GdgBMr6zJzSURcCezexW67U4zYqnY5RTvdWzsBbVQluDLzvoh4rDz+vzqJdQQwomrVGID29vY+j36u7N/TcYZsvz1Df/ELltx8M4sH+Yjr3taZlrHOamed1c46q10968x6b30msFrH7rz8TuvlFCOxemsSMIHlO7QzI+LG8vidJrCaoVO7ww47cPXVV3PLLbdw6KGH9umcA50NZ+2ss9pZZ7Wzzmo3CDq1awJDKW4UVXsa2LKLfSZ0UX5CDeedACzMzBdrOM4XgK91XDl58mRGjRpVw6m7dsUVV3S7ffX583kdsOCGG5h8ySV1OedA11Od6eWss9pZZ7WzzmpXjzqbO3duHSJRMzOB1Tq66tCOjYiRmTmvl8eo7NfxON11jBveqa2MKLvqqqu4xE4tYMO5Iqyz2llntbPOamentimcyPIjv8YAT+y7776MHTu2Twdub2/niiuuYJ999qGtra3rgq9/PfmlLzHy+efZf8cdYUItObvW0us601LWWe2ss9pZZ7WrZ53NmjWrTlGpWZnAUj00vFO76aab8r3vfY/HHnuMfffdl2HDBu9/2jactbPOamed1c46q90g6NTOABYDHd/+tw4wvYt9ptdYvqtjDI+I1TqMwuryOJm5AFhQ+V25cdTW1la3/557PNb48bDllnDvvbTdeSdssEFdzjuQ1bP+BwvrrHbWWe2ss9rVo86s89Y3eP+V33q66tDO6uXoq8oxKvtN63Cc27vaqRk6tVtttRVjxoxh9uzZPPTQQ2yzzTZ1Oe9AZsNZO+usdtZZ7ayz2rVqpzYzF0bErcDeFG8TJCKGlL9P72K3G8rtp1St26dc31u3Au3lcX5XnncLYMMaj7Py7bRT8RbCW2+F/fdvdDSSJGkl8i2EraPSoa1Wa4d2KkUSa+lxImIsxdsIm7pDO2TIEHYo3050m28nkiQNHCcDH42ID0XEVsCPKd4afDZARPwiIk6sKn8qsF9EHBsRW0bE14GdqUp4RcTqEbE98Mpy1RYRsX1ETIBifkvgLODkiHhjROxUnu+GzHzZBO5NZccdi89//7uxcUiSpJXOBFaTiohVy87m9uWqSeXvDcvtJ0bEL6p2+QmwcUScVHZoPwm8B/h+b4+ZmUlxR/fLEXFARGwD/AJ4ivLOcDOrJLD+badWkjRAZOYFwHHACRSjnbcH9svMynyUGwITq8pfDxwCfAy4AzgIODAz76o67AHAv4G/lr9/U/4+sqrMfwN/oRiB9Q+KG1jvqt+V9ZOyrcebVZIkDTo+Qti8dgaurvpdmWPqXOAwis7shpWNmTk1It5KkbD6NPAE8JHMvLyGYwKcRHHn9wxgNeA6io70/L5eUH+rJLBuv/32xgYiSVINMvN0unhkMDP37GTdRcBF3RzvHOCcHs45HziqXAaO7bcvPh99FJ5/HlZfvaHhSJKklccEVpPKzGuA6Gb7YV3ss8OKHrMsk8BXy2VA2b7s1N5+++1k5tK5uCRJUotYbTWYNAmmToXbb4e99mp0RJIkaSXxEUK1jFe+8pUMHz6cmTNnMnXq1EaHI0mS+kPlMUKnDJAkaVAxgaWW0dbWxtZbbw04D5YkSS2r8hihbb0kSYOKCSy1lOrHCCVJUgtyBJYkSYOSCSy1lMpE7nfccUeDI5EkSf2iMgLr/vth4cKGhiJJklYeE1hqKa985SsBuO+++xociSRJ6hfrrQerrgqLF8PDDzc6GkmStJKYwFJL2WKLLQCYMmUKC70rK0lS64mAsr3HG1aSJA0aJrDUUtZdd11WXXVVFi9ezMPelZUkqTVVElj339/YOCRJ0kpjAkstJSKWjsK6306tJEmtacsti09HYEmSNGiYwFLLMYElSVKLcwSWJEmDjgkstZwty7uyTuQuSVKLqh6BldnYWCRJ0kphAkstZ/PNNwfgwQcfbHAkkiSpX2y6afH54ovFIkmSWp4JLLWcSZMmAfDII480NhBJktQ/Ro2Ctdcuvk+d2thYJEnSSmECSy3nFa94BQBPPfUUCxYsaGwwkiSpf5Q3rExgSZI0OJjAUstZa621GDVqFJnJY4891uhwJElSf6gksBxxLUnSoGACSy0nIpaOwvIxQkmSWlTZ1jsCS5KkwcEEllpSZR6sqXZqJUlqTT5CKEnSoGICSy3JEViSJLU4HyGUJGlQMYGlluQILEmSWlx1AiuzoaFIkqT+ZwJLLamSwHIEliRJLWqDDSAC5s6FZ59tdDSSJKmfmcBSS6o8QugILEmSWtSIEbDeesV3b1hJktTyTGCpJW244YYAPP300yxYsKDB0UiSpH6xwQbF5+OPNzYOSZLU70xgqSWtscYarLLKKgA89dRTDY5GkiT1CxNYkiQNGiaw1JIigvXXXx+Ax+3USpLUmkxgSZI0aJjAUsvaoOzUmsCSJKlFmcCSJGnQMIGllmUCS5KkFleOtuaJJxobhyRJ6ncmsNSyTGBJktTiHIElSdKgYQJLLcsEliRJLa6SwHrqKVi0qLGxSJKkfmUCSy3LSdwlSWpx66wDbW2wZAlMm9boaCRJUj8ygaWWVRmB9YTzYkiS1JqGDIH11iu+e8NKkqSWZgJLLauSwJoxYwbz5s1rcDSSJKlfOA+WJEmDggkstazVVluN0aNHA47CkiSpZZnAkiRpUDCBpZYVEU7kLklqehFxVEQ8EhHzI+LGiNi1h/IHR8R9Zfk7I2L/DtsjIk6IiGkRMS8iroyIzTqU2Twi/hgRMyJiVkRcFxFv7I/r63cmsCRJGhRMYDWpiHh9RPw5Ip6KiIyIA3uxz54RcVtELIiIhyLisE7KdNtJjogJEfHLiJgeEXPK4727fle2cjmRuySpmUXEe4GTgW8AOwJ3AJdHxNpdlN8DOB84C9gBuBi4OCK2rip2PPAp4EhgN2BOecxVqsr8BRgG7AXsVJ73LxExoW4Xt7KYwJIkaVAwgdW8RlN0Jo/qTeGImAT8Fbga2B44BfhZRLy5qkxvOsm/ALYADgC2AX4PXBgRO/TtchrDEViSpCb3GeDMzDw7M++hSDrNBY7oovyngcsy87uZeW9mfgW4DTgaitFXwDHA/2TmHzPzP8ChwLrAgWWZNYHNgG9n5n8y80Hg88AoYGsGGhNYkiQNCiawmlRmXpqZX87MP/RylyOBqZl5bNmhPR34LfDfVWV600neAzgtM2/KzCmZ+T/AixR3ZwccE1iSpGYVEcMp2tcrK+syc0n5e/cudtu9unzp8qryk4AJHY45E7ixqsxzwP3AoRExOiKGAR8HngFu7cMlNUYlgeV8l5IktbRhjQ5AddNVh/YUWK6TfGJlY2YuiYiOneTrgfdGxF8pElfvAVYBrumnuPuVCSxJUhNbExgKPN1h/dPAll3sM6GL8hOqttNdmczMiHgTxeOHs4ElFMmr/TLzhc5OGhEjgBFVq8YAtLe3097e3kWovVPZf4WPM2ECbUA+/TSL5syB4cP7FM9A0Oc6G4Sss9pZZ7WzzmpXzzqz3lufCazW0VWHdmxEjATG07tO8nuACyjuzi6iGKH1zsx8qKsTN3OnduLEiUCRwBosf6HZcNbOOquddVY766x2dmr7R/mY4Q8pklavA+YBHwH+HBG7ZOa0Tnb7AvC1jisnT57MqFGj6hLXFVdcsWI7ZvK24cMZunAh1/z618xdZ526xDMQrHCdDWLWWe2ss9pZZ7WrR53NnTu3DpGomZnAUkffBFYD3gTMoJgv48KIeF1m3tnFPk3bqa2MvJo6dSqXXHJJXWIZKGw4a2ed1c46q511VrsW7tTOABYDHTMu6wDTu9hneg/lp1etm9ahzO3l972AtwHjM3NWue6TEbEP8CHg252c90SKeTQrxgBP7LvvvowdO7aLUHunvb2dK664gn322Ye2trYVOsaQDTeEhx7ijZtuSr7udX2KZyCoR50NNtZZ7ayz2llntatnnc2aNavnQhrQTGC1jq46tLMyc15ELKaHTnJEbEIxCezWmXl3uf2OiHgdxWTyR3Zx7qbt1M6aNYv/+q//Ys6cObzuda9jzJgxfYpnILDhrJ11VjvrrHbWWe1avVObmQsj4lZgb4rH+YiIIeXv07vY7YZy+ylV6/Yp1wNMpWjX96ZMWEXEWIq3Ef64LFO5u7Skw7GX0MX8qJm5AFhQ+V0M4oK2tra6/ffcp2OVCaxhTz0Fg+j/r3rW/2BhndXOOquddVa7etSZdd76TGC1jhuA/TusW9qh7WUnuasO7WK6mfC/mTu1a6yxBuPGjWPmzJlMnz6d1VdfvS7xDAQ2nLWzzmpnndXOOqtdi3dqTwbOjYhbgJso3iA4GjgbICJ+ATyZmV8oy58K/D0ijqV4+/D7gJ2Bj8HS+a1OAb4cEQ9SJLS+CTxF2f5T9A1eKM97AsUjhB+lmAD+r/14rf1nww2LT+e8lCSpZfkWwiYVEatGxPYRsX25alL5e8Ny+4llp7biJ8DGEXFSRGwZEZ+kmM/q+1VlTgY+GhEfioitKO7ELu0kA/cBDwE/jYhdI2KTsoO8D8s6vQPOhmWn9rHHHmtwJJIkLS8zLwCOA06gGDG1PcVk6pU5KzcEJlaVvx44hCJhdQdwEHBgZt5VddiTgNOAM4CbgVXLY84vjzED2K9cfxVwC/Ba4B2ZeUd/XGe/qySwHn20sXFIkqR+4wis5rUzcHXV78ojeucCh1F0ZjesbMzMqRHxVoqE1aeBJ4CPZOblVWUuiIi1KDrJEyg6yks7yZnZHhH7U8x98WeKju1DwIcyc8BOILXRRhtx5513msCSJDWlzDydLh4ZzMw9O1l3EXBRN8dL4Kvl0lWZW4A31xpr06oksGzrJUlqWSawmlRmXgNEN9sP62KfHXo4bped5HL7g8C7exnmgFAZgfWod2UlSWpNG21UfJrAkiSpZfkIoVqejxBKktTiqh8hzGxsLJIkqV+YwFLL26i8K2sCS5KkFrXBBsXnSy/Biy82NBRJktQ/TGCp5fkIoSRJLW7kSFh77eK7N6wkSWpJJrDU8ioJrCeeeILFixc3OBpJktQvfBOhJEktzQSWWt7EiRMZNmwYixcvZtq0aY0OR5Ik9QffRChJUkszgaWWN3ToUNZff33AxwglSWpZJrAkSWppJrA0KPgmQkmSWlz50hYfIZQkqTWZwNKgYAJLkqQW5wgsSZJamgksDQoblXdlfYRQkqQWZQJLkqSWZgJLg4IjsCRJanGVBNa0abBwYWNjkSRJdWcCS4OCCSxJklrcWmvBKqtAJjzxRKOjkSRJdWYCS4NCJYHlI4SSJLWoCB8jlCSphZnA0qBQSWDNmjWLmTNnNjgaSZLULyoJLG9YSZLUckxgaVBYddVVWX311QFHYUmS1LLKl7bwyCMNDUOSJNWfCaw6iYi2iNggIraIiNUbHY9ebtKkSQBMmTKlwZFIkgYi2/oBYOONi8+HH25sHJIkqe5MYPVBRIyJiE9ExN+BWcAjwL3AsxHxaEScGRG7NDRILbXJJpsA8LCdWklSL9nWDzBlW28CS5Kk1mMCawVFxGcoOrGHA1cCBwLbA5sDuwPfAIYBkyPisojYrCGBaikTWJKkWtjWD0AmsCRJalnDGh3AALYL8PrMvLuL7TcBP4+IIyk6vq8DHlxZwenlTGBJkmpkWz/QVBJYTz8NL70Eq67a2HgkSVLdmMBaQZn5foCIGAq8HfhbZs7upNwC4CcrOTx1wgSWJKkWtvUD0PjxsPrq8PzzMGUKbLttoyOSJEl14iOEfZSZi4HzgbUaHYu6V0lgPfrooyxatKjB0UiSBgrb+gHGxwglSWpJJrDq42ZgUqODUPfWW289RowYwaJFi3jssccaHY4kaWCxrR8oTGBJktSSTGDVx2nAtyJig0YHoq4NGTKESZOKf3v4GKEkqUa29QOFCSxJklqSCaz6uIBiote7I+JXEfGRiNgpIoY3OjAtz3mwJEkryLZ+oDCBJUlSS3IS9/qYBGxH8Wrt7YAvAK8AFkXE/ZnpDKJNwgSWJGkF2dYPFCawJElqSSaw6iAzHwUeBf5UWRcRYyg6uXZom4gJLEnSirCtH0AqCaxHH4X2dmhra2w8kiSpLkxg9ZPyNdvXlouaRCWB9dBDDzU4EknSQGdb36QmToSRI2HePHjkEdhss0ZHJEmS6sA5sDSobLHFFgA8+OCDLFmypMHRSJKkuhsyBDbfvPh+332NjUWSJNWNCSwNKq94xSsYPnw48+fP59FHH210OJIkqT9stVXxaQJLkqSWYQJLg8qwYcPYvLwre5+dWkmSWtOWWxaf997b2DgkSVLdmMDqZxGxJCKuioidGh2LCluWndp77dRKkurAtr4JOQJLkqSWYwKr/x0B/AP4YaMDUWGrslPrCCxJUp3Y1jeb6hFYmY2NRZIk1YVvIexnmXlO+fXrDQxDVRyBJUmqJ9v6JrTZZhABL74IzzwD66zT6IgkSVIfOQKrSUXE6yPizxHxVERkRBzYi332jIjbImJBRDwUEYd1UuaoiHgkIuZHxI0RsWsnZXYvH4WYExGzIuIfETGyPlfWeJUEliOwJEnNoDdtc4fyB0fEfWX5OyNi/w7bIyJOiIhpETEvIq6MiM06Oc5by/PNi4gXIuLiOl9a44wcCZMmFd+9YSVJUkswgdVHEbFmRBwfEX+IiBvK5Q8R8dmIWKsPhx4N3AEc1cs4JgF/Ba4GtgdOAX4WEW+uKvNe4GTgG8CO5fEvj4i1q8rsDlwGTAZ2BXYBTgeW9OFamsoWW2wBwIwZM5gxY0aDo5EkNbt+bOt71TZ3KL8HcD5wFrADcDFwcURsXVXseOBTwJHAbsCc8pirVB3n3cAvgbOB7YDXAOf15VqaTuUxQm9YSZLUEkxg9UFE7AI8QNFJnEkx/8U/yu+fAu6LiJ1X5NiZeWlmfjkz/9DLXY4EpmbmsZl5b2aeDvwW+O+qMp8BzszMszPznnKfuRRzd1R8H/hBZn47M+/OzPsz88LMXLAi19GMRo8ezYYbbgg4CkuS1L3+bOtLvWmbq30auCwzv1u2918BbgOOLuMN4BjgfzLzj5n5H+BQYF3gwLLMMOBU4LOZ+ZPMfCAz78nMC/twHc2nMpG7I7AkSWoJzoHVN6cBFwFHZi4/Q2jZgfxJWWb3lRDL7sCVHdZdTjESi4gYDuwEnFjZmJlLIuLKSnzl3d7dgF9HxPXAJsB9wJcy87r+voCVaauttuKxxx7j3nvv5bWvfW2jw5EkNa9+a+t70zZ3YneKEVvVLqdMTgGTgAlU9Qkyc2ZE3Fju+xuKkV7rAUsi4t9l+dspElp3dRHrCGBE1aoxAO3t7bS3t/d0qd2q7N/X43QUm23GMGDJ3XezuM7HbrT+qrNWZp3VzjqrnXVWu3rWmfXe+kxg9c12wGEdO7QAmZkR8X3g3ysplgnA0x3WPQ2MLeevGg8M7aJMOcaejcvPrwPHUXRmDwX+FhFbZ+aDnZ14IHZqt9xySy6//HLuuOOOlvuLzoazdtZZ7ayz2llntWuSTm1/tvVr0nPb3FFX7f2Equ30UKa6vf8M8AhwLHBNRGyemc93ct4vAF/ruHLy5MmMGjWqi1Brc8UVV9TlOBXjZ83i9cDCW2/l8ksuqeuxm0W962wwsM5qZ53VzjqrXT3qbO7cuXWIRM3MBFbfTKeYJ6qr59B25eUdyGZWeaT0p5l5dvn93xGxN8WjDF/oYr8B16mt/Dvkmmuu4RI7tSpZZ7WzzmpnndWuwZ3aVmvrYVl7/7+Z+TuAiDgceAI4GPhpJ/ucyPIjv8YAT+y7776MHTu2T8G0t7dzxRVXsM8++9DW1tanYy3nDW8gP/c5VnnxRfbfaaeWehNhv9VZC7POamed1c46q10962zWrFl1ikrNygRW3/wfcEZE7AT8jWUd2HWAvYGPUoxkWhmml+ettg4wKzPnRcRiYHEXZaaX36eVn/d0KHMvsGE35x5wndp1112XH/zgBzz55JO85S1voXgKpDXYcNbOOquddVY766x2TdKp7c+2fgY9t80dddXeT6/aXlk3rUOZ28vvL2vvM3NBREyhi/a+nAtz6XyYlXazra2tbv891/NYAKy2Gmy2GTzwAG333APrr1+/YzeJutfZIGCd1c46q511Vrt61Jl13vpMYPVBZv4wImZQTJT+SYrHAKDojN5K8cjBypoQ9QZg/w7r9inXk5kLI+JWis72xQARMaT8fXpZ/hHgKWCLDsfZHLi0qxMPxE7ttttuy7Bhw3jxxReZPn360kndW4kNZ+2ss9pZZ7WzzmrXyE5tf7b1vWybO7qh3H5K1bql7T0wlSKJtTdlwioixlLMcfnjssytFO32FsB1ZZk24BXAoytyLU1r223hgQfgjjtg330bHY0kSeoDE1h9lJkXABeUHb81y9UzMrNPE3ZExKrAplWrJkXE9sDzmflYRJwIrJeZh5bbfwIcHREnAT8H9gLeA7y16hgnA+dGxC3ATRRvKRpN8Qrtylwe3wW+ERF3UHR8P0QxD8dBfbmeZjNixAi22mor7rzzTu64446WTGBJkuqjv9r6Urdtc0T8AngyMyuP8Z8K/D0ijgX+CrwP2Bn4WBlrRsQpwJcj4kGKhNY3KW5QXVyWmRURP6Fo7x+nSFp9tjz+RXW4puax3Xbw298WCSxJkjSgmcBaQRGxYWY+VvlddmKndVN+vcx8soZT7AxcXfW78ojeucBhwESqhvln5tSIeCvwfYpXbD8BfCQzL68qc0FErAWcwLI3Du2XmU9XlTklIlYpj7M6cAewT2Y+XEPsA8J22223NIH19re/vdHhSJKazEpo63vTNm8ILKkqf31EHAL8D/At4EHgwA5vDzyJIgl2BrAaxSir/TJzflWZzwKLgF8CI4Ebgb0y84Va4m96221XfP7nP42NQ5Ik9dmQnouoCzdHxE8jYpeuCkTEuIj4aETcBby7loNn5jWZGZ0sh5XbD8vMPTvZZ4fMHJGZm2TmOZ0c9/TM3Kgss1tm3thJmW9n5gaZOToz98jM62qJfaDYruzU/sdOrSSpc/3a1ld01zZn5p6Vtr9q3UWZuUVZfuvMvKTD9szMr2bmhMxcJTPflJkPdCjTnpnHZeY6mTk2M/fJzLtXJP6mVklg3XsvLFjQfVlJktTUHIG14l4JfAm4IiLmU8wn8RQwHxhfbn8VcBtwfMfOpRpv2223BeAOHyuQJHXOtn6g22CDYjL3F18skljbb9/ggCRJ0opyBNYKysznMvMzFI/yHUUxhH9NYLOyyK+BnTJzdzu0zakyAuvBBx9kzpw5DY5GktRsOrT1R2NbP/BEFBO5A9x+e0NDkSRJfeMIrD7KzHkRMRv4VmY+0+h41HvrrLMOEydOZNq0adx+++285jWvaXRIkqQmlJnzgN+WiwaanXeGf/wDbrkFDjus0dFIkqQV5Ais+rgYmBYRT0XEJRHxvxGxZ2NDUm/suuuuANx0000NjkSS1OwiYtuIOKqc8+qVjY5HvVS29djWS5I0oJnAqo8xwHYUb/S5C9gNuCwiro6I0Q2NTN0ygSVJ6o2I+DTFGwL/F/g2cFdE3BER2zcyLvXCLuUc/Lff7kTukiQNYCaw6mONzLwrM3+dmcdn5psoXnvdBnylwbGpGyawJEldiYgjImLHiBhBMZn754HxmbkGsDFwKXBtROzRyDjVg0mTYI01oL0dfHGLJEkDlgms+phePj54aUR8OyLeD6wBHAMc0djQ1J2dd94ZgClTpjBjxowGRyNJajLHATcCL1G067sAn46INwAvZObnKZJa/9e4ENWjiGWPEd58c2NjkSRJK8wEVn1sRvF2ohuBLSgeL7gbuA5YIyJ+FRH/FRGva2CM6sRqq63G5ptvDsDNdmolSVUy85UU0wTsAbQDS4D3AZcAz0fEFOCdwE4R8daIeEWjYlUPnAdLkqQBzwRWHWTmw5n5+8z8ema+MzM3BlYD3gtEuXwIuLyBYaoLlccITWBJkjrKzPmZeTPwT+COzHw1RVJrG+DLwEMUUwb8ApgSEbMaFqy6ZgJLkqQBb1ijA2gFEfE8xcSud5TLnRSPG+wPPJyZHyjLDW1UjOrarrvuyq9+9StuvPHGRociSWpexwLXRMTGwE8o2vvHgR2BpzJz/YhYH9i6gTGqK5WJ3O+/H158EVZbrZHRSJKkFWACqz6OoHgL4XbAO4BXlOvnAu+pFMrMxSs9MvVot912A+CGG25gyZIlDBniwERJ0vIy8/aI2IkiefUvitHVAIso57vMzCeAJxoTobq11lqw6abw0ENw/fWw//6NjkiSJNXIBFYdZObFwMWV3xExBpgIPJmZcxoUlnpphx12YNSoUbzwwgvcfffdbLPNNo0OSZLUhDLzYWCfiFgHeDUwHLihTFyp2b3+9UUC69prTWBJkjQAOdSkH2Tm7Mx8wOTVwNDW1sYeexRvQL/22msbHI0kqdll5tOZ+cfMvMjk1QDyuvJdOv/4R2PjkCRJK8QElgS8/vWvB+AfdmolSWpNZVvPzTfDvHmNjUWSJNXMBJYEvK68K3vttdeSmQ2ORpIk1d2kSbDeetDeDr64RZKkAccElkQxkXtbWxtPPfUUU6ZMaXQ4kiSp3iJ8jFCSpAHMBJYEjBw5kl133RVwHixJklpW5TFCE1iSJA04JrCkUmUerKuvvrrBkUiSpH7xhjcUn//8p/NgSZI0wJjAkkpvetObAJg8ebLzYEmS1Iq22qqYB2v+fLjuukZHI0mSamACSyq95jWvYdSoUUyfPp277rqr0eFIkqR6i4B99y2+T57c2FgkSVJNTGBJpREjRrDnnnsCcPnllzc2GEmS1D8qCSzbekmSBhQTWFKVfctO7WTvykqS1Jre9KZiJNadd8K0aY2ORpIk9ZIJLKlKJYH1j3/8g3lO7ipJUutZc03Yaafi+xVXNDYWSZLUayawpCpbbrkl66+/PgsWLOAfvmJbkqTWVHmM8LLLGhuHJEnqNRNYUpWIYL/99gPgL3/5S4OjkSRJ/WL//YvPSy6B9vbGxiJJknrFBJbUwTve8Q4ALr74YjKzwdFIkqS6e/WrYe21YeZM+PvfGx2NJEnqBRNYUgd77703o0eP5oknnuC2225rdDiSJKnehg6FAw4ovv/xj42NRZIk9YoJLKmDkSNHLn2M8OKLL25sMJIkqX+UI665+GJwxLUkSU3PBJbUiQMPPBAwgSVJUsvae28YPRqeeAIccS1JUtMzgSV14q1vfStDhw7lrrvu4qGHHmp0OJIkqd5GjoRyxDV/+ENjY5EkST0ygSV1Yvz48bzxjW8E4KKLLmpwNJIkqV+8+93F529+42OEkiQ1ORNYUhfe//73A/DrX//atxFKktSKDjgARo2Chx+Gm25qdDSSJKkbJrCkLrzrXe9i+PDh3H333fznP/9pdDiSpBYVEUdFxCMRMT8iboyIXXsof3BE3FeWvzMi9u+wPSLihIiYFhHzIuLKiNisi2ONiIjbIyIjYvs6XtbAMHo0lPNect55DQ1FkiR1zwRWk4qI10fEnyPiqbJTeWAv9tkzIm6LiAUR8VBEHNZJmV51ksvO76W9PXcrWm211Xjb294GwHl2aiVJ/SAi3gucDHwD2BG4A7g8ItbuovwewPnAWcAOwMXAxRGxdVWx44FPAUcCuwFzymOu0skhTwKeqsvFDFQf+EDx+ZvfwKJFjY1FkiR1yQRW8xpN0Yk9qjeFI2IS8FfgamB74BTgZxHx5qoytXSSjwEG/XNzHyg7teeffz5LlixpcDSSpBb0GeDMzDw7M++hSDrNBY7oovyngcsy87uZeW9mfgW4DTgaihtQFG34/2TmHzPzP8ChwLrAgdUHioi3APsCx9X9qgaSffaBNdeEZ56Bq65qdDSSJKkLJrCaVGZemplfzszevhbnSGBqZh5bdmhPB34L/HdVmV51kstHCI7tuH4w2n///Rk3bhyPP/44f//73xsdjiSphUTEcGAn4MrKusxcUv7evYvddq8uX7q8qvwkYEKHY84Ebqw+ZkSsA5wJfJCiLzB4tbXBe95TfD/nnIaGIkmSujas0QGobrrq0J4Cy3WST6xszMwlEbFcJzkiRgHnAUdl5vTiRm73ImIEMKJq1RiA9vZ22tvbV+Ralqrs39fjrKihQ4dy8MEH87Of/YwzzjiD1772tQ2JoxaNrrOByDqrnXVWO+usdvWssyat9zWBocDTHdY/DWzZxT4Tuig/oWo73ZUpR2mdA/wkM2+JiFf0FGgrt/UAHHoobT/6Efm737Fo+nRYY43GxdILTVFnA4x1VjvrrHbWWe0GQVuvOjKB1Tq66tCOjYiRwHh610n+PnB9Zv6xhnN/Afhax5WTJ09m1KhRNRyma1dccUVdjrMittyyqJ7f/va37L///owdO7ZhsdSikXU2UFlntbPOamed1a4edTZ37uAeZNTBf1EkoE7sqWCVlm7rAd6w8casNmUK933pS0w54ICGxtJbja6zgcg6q511VjvrrHa29eoNE1haKiIOAPaimBS2FidSzK1VMQZ4Yt999+1zsqe9vZ0rrriCffbZh7a2tj4dqy9+/etf8+9//5vp06fzvve9r2Fx9Eaz1NlAYp3VzjqrnXVWu3rW2axZs+oUVV3NABYD63RYvw4wvYt9pvdQfnrVumkdytxeft+LYvT1gg4jrW+JiF9n5oc6OW/Lt/VDnnwSjjqKrf/5T7b88Y+hF6PQG6VZ6mwgsc5qZ53Vzjqr3SBo61VHJrBaR1cd2lmZOS8iFtNzJ3kvYBPgxQ4d2t9FxLWZuWdnJ87MBcCCyu/Kvm1tbXX7i7uex1oRRx55JB//+Mc566yzOO644+jNo5WN1ug6G4iss9pZZ7WzzmpXjzprxjrPzIURcSuwN8XbBImIIeXv07vY7YZy+ylV6/Yp1wNMpWjX96ZMWEXEWIq3Ef64LPMp4MtV+69LMe3Aeynmyuos1pZv6/l//w+OP564/37abrwRXve6xsXSSw2vswHIOquddVY766x2rdrWq76cxL11VDq01ZZ2aDNzIXBrdZmqTnKl0/ttYFuKtxhWFigmgj+8X6IeIN7//vczevRo7r//fq655ppGhyNJah0nAx+NiA9FxFYUSabRwNkAEfGLiKh+1O9UYL+IODYitoyIrwM7Uya8MjMpkltfjogDImIb4BfAU5RJssx8LDPvqizAA+WxH87MJ/r3cpvY2LHw/vcX33/4w8bGIkmSXsYEVpOKiFUjYvvyjYAAk8rfG5bbT4yIX1Tt8hNg44g4qezQfhJ4D8WcVhXddpIzc3p1h7bs1AI8lplT++9qm9+YMWP44Ac/CMD3v//9HkpLktQ7mXkBcBxwAsWIqe2B/TKzMmflhsDEqvLXA4cAHwPuAA4CDqxqswFOAk4DzgBuBlYtjzm/P6+lJRx9dPH529/CY481NhZJkrQcE1jNa2fg3+UCRfLp3xQdXCg6sxtWCpcJprdSjLq6AzgW+EhmXl5VpqdOsrpxzDHHAPDnP/+Z+++/v7HBSJJaRmaenpkbZeaIzNwtM2+s2rZnZh7WofxFmblFWX7rzLykw/bMzK9m5oTMXCUz35SZD9CFzHwkMyMzb6/3tQ04220He+8NixfDqac2OhpJklTFBFaTysxrys5kx+WwcvthHeekKvfZoezQbpKZ53Ry3C47yV3EEZl5cf2ubODaYostePvb3w44CkuSpJb1mc8Un2eeCTNnNjYWSZK0lAksqQbHHnssAOeeey7PPvtsg6ORJEl1t99+sNVWMHs2/OxnjY5GkiSVTGBJNXj961/PTjvtxPz58znVRwskSWo9Q4YsG4X1ve/BvHmNjUeSJAEmsKSaRARf/OIXAfjBD37Ac8891+CIJElS3R16KGy4IUybBmec0ehoJEkSJrCkmh144IFst912zJ49m+9973uNDkeSJNXb8OHwpS8V37/9bUdhSZLUBExgSTUaMmQI3/jGNwA47bTTmDFjRoMjkiRJdXfYYbDRRjB9OvzkJ42ORpKkQc8ElrQCDjjgAHbccUdeeuklvvOd7zQ6HEmSVG/Dh8OXv1x8P/FEmDWrsfFIkjTImcCSVkBE8M1vfhMo5sKaMmVKgyOSJEl196EPweabw7PPFkksSZLUMCawpBX0lre8hTe96U0sXLiQz33uc40OR5Ik1VtbG/zf/xXfv/99eOSRhoYjSdJgZgJLWkERwfe+9z2GDBnCb3/7W6699tpGhyRJkurtbW+DvfaCBQvg859vdDSSJA1aJrCkPth222358Ic/DMAxxxzD4sWLGxyRJEmqqwg4+eTi84IL4O9/b3REkiQNSiawpD765je/ybhx47jttts4/fTTGx2OJEmqt+22g49/vPj+8Y/D/PmNjUeSpEHIBJbUR+uss87SNxF+6Utf4rHHHmtwRJIkqe5OPBEmTID773dCd0mSGsAEllQHH/3oR3nNa17DnDlzOOqoo8jMRockSZLqabXV4LTTiu8nngj33NPQcCRJGmxMYEl1MGTIEM444wza2tr4y1/+wnnnndfokCRJUr29+93w9rdDezscdljxKUmSVgoTWFKdvPKVr+SrX/0qAJ/85Cd59NFHGxyRJEmqqwj40Y9g/Hi4+WY44YRGRyRJ0qBhAkuqo89//vPssccezJo1iw9+8IO+lVCSpFaz/vrw058W37/1LbjuusbGI0nSIGECS6qjYcOG8ctf/pIxY8Zw7bXX8u1vf7vRIUmSpHo7+GA49FBYsgQ++EF44YVGRyRJUsszgSXV2cYbb8xp5SSvX/3qV7nyyisbHJEkSaq7006DSZPgkUfg//2/IpklSZL6jQksqR8ceuihHH744SxZsoT3ve99zoclSVKrGTsWfvc7WGUVuOQS58OSJKmfmcCS+kFE8MMf/pCddtqJ5557jne9613Mmzev0WFJkqR62mGHZfNhfeMb8Je/NDYeSZJamAksqZ+MHDmS3/3ud6yxxhrcdtttHHbYYSzx8QJJklrLoYfCUUcV39/3PrjttsbGI0lSizKBJfWjjTbaiN/+9re0tbVx4YUXcvzxxzc6JEmSVG8nnwxvehPMmQNvfSs4dYAkSXVnAkvqZ3vuuSdnn302AN/73vf4wQ9+0OCIJElSXQ0fDr/9LWyzDUyfDm95i28mlCSpzkxgSSvBBz7wAb71rW8BcMwxx/CLX/yiwRFJkqS6GjeumMx9vfXg3nvhzW+GmTMbHZUkSS3DBJa0knz+85/n6KOPJjM5/PDDOe+88xodkiRJqqf114fLLoM11oCbby5GYs2e3eioJElqCSawpJUkIjj11FP56Ec/ypIlS/jgBz/IhRde2OiwJElSPW29NVx5JYwfDzfcUMyJ9dJLjY5KkqQBzwSWtBINGTKEn/zkJxx++OEsWbKEQw45ZOn8WJIkqUVsvz1ccUXxWOG118Jee8GMGY2OSpKkAc0ElrSSDRkyhDPPPJPDDz+cxYsXc8QRR3DSSSc1OixJklRPO+1UJLEqjxO+7nXw2GONjkqSpAHLBJbUAEOHDuWss87is5/9LACf+9znOO6441iyZEmDI5MkSXWzyy7FCKz114f77oPXvAb+859GRyVJ0oBkAktqkIjgpJNO4rvf/S4A3/ve9zjwwAOZNWtWgyOTJEl1s9VWcP31sOWW8MQTsMce8Mc/NjoqSZIGHBNYUoMdd9xx/PrXv2bEiBH8+c9/Zo899uDhhx9udFiSJKleNtgA/vlP2HtvmDMHDjwQvvUtyGx0ZJIkDRgmsKQmcMghh/CPf/yDiRMncvfdd7Prrrvy17/+tdFhSZKkell9dbj0Ujj66OL3l74E7343vPBCY+OSJGmAMIHVpCLi9RHx54h4KiIyIg7sxT57RsRtEbEgIh6KiMM6KXNURDwSEfMj4saI2LVq2+oRcVpE3B8R8yLisYj4QUSMq+/VqTO77rort9xyC7vssgvPP/88b3vb2zj22GNZuHBho0OTJPWj7trmLsofHBH3leXvjIj9O2yPiDghIqaV7fmVEbFZ1fZXRMRZETG13P5wRHwjIob31zWq1NYGp50GP/1p8f0Pf4AddoAbbmh0ZJIkNT0TWM1rNHAHcFRvCkfEJOCvwNXA9sApwM8i4s1VZd4LnAx8A9ixPP7lEbF2WWTdcjkO2Bo4DNgPOKuvF6PeWXfddbn22mv51Kc+BcDJJ5/Ma17zGh566KEGRyZJ6g+9aJs7lt8DOJ+ibd4BuBi4OCK2rip2PPAp4EhgN2BOecxVyu1bUvQBPw68Cvjvsuy36nlt6sbHPlbMi7XJJvDoo8UbCr/9bVi8uNGRSZLUtExgNanMvDQzv5yZf+jlLkcCUzPz2My8NzNPB35L0Smt+AxwZmaenZn3lPvMBY4oz3lXZr47M/+cmQ9n5lXAl4C3R8Swul2cujVixAhOPfVULr74YsaPH88tt9zCdtttx6mnnupbCiWp9XTbNnfi08Blmfndsr3/CnAbcDQUo6+AY4D/ycw/ZuZ/gEMpblAdCJCZl2Xm4Zk5OTOnZOafgP8D3tVvV6mX23lnuO02eN/7isTVF75QvKXw3nsbHZkkSU3JpETr2B24ssO6yylGYlE+FrATcGJlY2YuiYgry327Mg6YlZmLuioQESOAEVWrxgC0t7fT3t5ewyW8XGX/vh5nINp///255ZZb+PCHP8w111zDMcccw4UXXsgZZ5zB5ptv3uV+g7nOVpR1VjvrrHbWWe3qWWfNWO8r2DbvTjFiq9rllMkpYBIwgao+QWbOjIgby31/08VxxwHPdxOrbX1/GDkSzj2XeNObGPqZzxA33kjusANLvvIVlnzmMzCs8676oK6zFWSd1c46q511VrtWb+tVX5G+/aTpRUQC78zMi7sp8wBwdmaeWLVuf4rHCkcB44EngT0y84aqMicBb8jM3To55prArcCvMvNL3Zz768DXOq4/77zzGDVqVI/Xp+4tWbKEyZMnc8455zB//nyGDx/Ou971Lt75zncyYsSIng8gSWLu3LkccsghAOMyc1aj4wGIiHWpvW1eCHwoM8+vWvdJ4GuZuU75iOE/gXUzc1pVmQuBzMz3dnLMTSna++My88wuYv06tvX9apUZM9juxz9mwq23AjDzFa/gPx/9KM+/6lUNjkySBoZmbOtVX47AUqciYixF8use4Os9FD+R5e8GjwGe2HfffRk7dmyf4mhvb+eKK65gn332oa2trU/HGsgqE7p/4hOf4Morr+Q3v/kN//rXv/jud7/LAQccQPHESME6q511VjvrrHbWWe3qWWezZtmP7UxErAdcBlzUVfKqZFu/Mnzwgyz65S8Z+tnPMu6RR3jdl77EkkMOYfGJJ8LEiUuLWWe1s85qZ53VzjqrnW29amECq3VMB9bpsG4disf/5kXEYmBxF2WmV6+IiDEUndnZFCO/uh2LmZkLgAVV+wPQ1tZWt7+463msgWrTTTdl8uTJXHjhhRx33HE88sgjHHzwwey777585zvfYfvtt1+uvHVWO+usdtZZ7ayz2tWjzpq0zmfQy7a5Slft/fSq7ZV10zqUub16p3IE2NXA9cDHugvUtn4l+vCH4R3vgC99Cc48kyHnnceQP/0JPv95OOYYGD16aVHrrHbWWe2ss9pZZ7Vr4bZedeQk7q3jBmDvDuv2KdeTmQspHg9YWiYihpS/qx9bGAtMBhYCB2Tm/P4NW7WICN773vdy33338cUvfpHhw4czefJkdthhBw455BDfVihJA0hv2+YOum3vgakUSazqY46leBthdXu/HnBNef7DM9O3hDSTNdeEn/4UbroJdtsNXnoJvvzl4q2FP/whLFzY6AglSVrpTGA1qYhYNSK2j4jty1WTyt8blttPjIhfVO3yE2DjiDgpIrYs58N4D/D9qjInAx+NiA9FxFbAj4HRwNnlMSvJq9HAh4GxETGhXIb24+WqRqNHj+Z///d/ufvuu3nf+94HwPnnn89WW23FUUcdxTPPPNPgCCVJvdRT2/yLiDixqvypwH4RcWzZ3n8d2Bk4HYpJrihe4PLliDggIrYBfgE8BVxcHrOSvHoMOA5Yq9Le9/O1qlY77wzXXw/nnQcbbwxPPw1HH82wbbZhw7/9zUSWJGlQMYHVvHYG/l0uUHRw/w2cUP6eCGxYKZyZU4G3UtyFvQM4FvhIZl5eVeYCio7qCRSPEWwP7JeZT5dFdqS4Q7sN8BDFoweVZYM6X5/qYNNNN+X888/ntttuY7/99mPRokWceeaZfOITn+CII47gnnvuaXSIkqRu9KJt3pCiza+Uvx44hOKRvzuAg4ADM/OuqsOeBJwGnAHcDKxaHrMyqnofYFOKUVpPsHx7r2YzZAi8//1w773wox/BhAnE1KnscNppDNtqK/jBD2Du3EZHKUlSvzOB1aQy85rMjE6Ww8rth2Xmnp3ss0NmjsjMTTLznE6Oe3pmblSW2S0zb+zFOSMzH+nfK1Zf7LDDDlx66aVcc801vPGNb2Tx4sX86le/4lWvehXvfOc7ueGGG/CNo5LUnHpom/estP1V6y7KzC3K8ltn5iUdtmdmfjUzJ2TmKpn5psx8oGr7OV219/1+sVpxw4fDJz4BDz3E4hNPZP748cTjj8OnPw0bbQTf/CY4AluS1MJMYEkt5A1veAOXX345J510Eu94xzsAuPjii9ljjz3YeeedOfvss5k3b16Do5QkSSts9GiWHHssV/z0pyz+4Q+LRwtnzICvfhU22AAOPbSYO0uSpBZjAktqQZtvvjkXXXQR99xzD4cffjgjRozgtttu44gjjmD99dfn+OOPd8J3SZIGsCXDh7Pkox+F++8v5sjabbdiTqxf/rL4vttucM45xQTwkiS1ABNYUgvbaqut+PnPf84TTzzBd77zHTbaaCOef/55vvvd77LZZpvxute9jrPOOotZs2Y1OlRJkrQihg0r5sj617+KkVeHHlo8bnjTTXD44TBhQvH597/DEl82KUkauExgSYPAmmuuyfHHH8/DDz/MH//4R/bbbz+GDBnCddddx0c+8hEmTJjABz/4QSZPnkx7e3ujw5UkSStil13g3HPhiSfgW9+CzTaDOXOKkVh77ln8/vrXiwnhJUkaYExgSYPI0KFDOeCAA7j00kt57LHH+Pa3v82WW27JvHnz+NWvfsWb3/xmJkyYwEc+8hEuu+wyk1mSJA1Ea60FX/hC8XjhddfBRz4CY8bAlCnwjW/AK18J22wDJ5wA993X6GglSeoVE1jSILXeeuvxuc99jnvuuYd//etffOITn2Cttdbi+eef56yzzuItb3kL66yzDkcccQR//OMfmTNnTqNDliRJtYiA17wGzjwTpk8v5sfaf39oa4O77oKvfQ222qpIZn3ta3DzzT5mKElqWiawpEEuIthtt9340Y9+xFNPPcXf/vY3PvGJT7D22mvzwgsvcPbZZ3PggQey+uqr8+Y3v5nTTjuNKVOmNDpsSZJUi1Gj4P/9P/jrX+Hpp4vHCquTWSecALvuCuuuC0ccAb/7HThHpiSpiZjAkrTUsGHD2GuvvZYms66++mqOPvpoJk2axMKFC5k8eTKf+tSn2GSTTdhqq6045phj+NOf/sTMmTMbHbokSeqt8ePhQx9alsw691w46KDiMcOnn4azzy5+r7km7LUX/O//wg03wKJFjY5ckjSIDWt0AJKa09ChQ9lzzz3Zc889+cEPfsB9993HX//6V/76179y3XXXcd9993Hfffdx6qmnMmTIEHbeeWf22msv9tprL17zmtcwatSoRl+CJEnqyfjxxZsLDz0UFi4s5sz6y1+K5cEH4eqriwWKBNcb3lAktfbeG7beGoZ4P1yStHKYwJLUo4hgq622YquttuK4445j5syZXHHFFfztb3/jb3/7Gw8++CA33XQTN910E9/+9rcZPnw4O+20E3vsscfSZcKECY2+DEmS1J3hw4vk1F57wcknwwMPwJVXwt/+ViSxXnhhWXILiuTX7rsX82ztsUfxCKI3sCRJ/cQElqSajRs3joMOOoiDDjoIgMcff5yrr756aULrySef5IYbbuCGG27ge9/7HgAbb7zx0mTWbrvtxtZbb83w4cMbeRmSJKk7m29eLJ/8ZDG5++23w1VXFQmtf/yjSGhdckmxAAwbBttvXyS0dt8ddtoJNtmkmExekqQ+MoElqc822GADDj30UA499FAyk4cffpjrr79+6XLXXXcxZcoUpkyZwq9+9SsAhg8fzjbbbMNOO+3EzjvvzE477WRSS5KkZjVkCOy4Y7Ecdxy0t8Mdd8A//wnXX198Pvkk3HJLsZx6arHfuHFFIquy7LwzbLyxSS1JUs1MYEmqq4hg0003ZdNNN+XQQw8FYObMmdx4441cf/31/POf/+SWW27hxRdf5NZbb+XWW2/ljDPOAJYltXbYYQe22WYbtt56a7bZZhvWWmutRl6SJEnqqK2tSEbtvDN8+tOQCY8/XiSy/vlPuOmmIsE1c2Yxauuqq5btu9pqsMMOsM02y5ZXvQpWXbVhlyNJan4msCT1u3HjxrHvvvuy7777ApCZTJ06dWkCq7K88MILS79XW3vttZcms7beemu23nprXvWqVzFmzJhGXI4kSeooAjbcsFje//5i3cKFcPfdcOuty5Y77oAXX1x+cviKSZNentTabDMYMWKlX44kqfmYwJK00kUEG2+8MRtvvDEHH3wwUCS1HnnkEW699VbuuOMO7rrrLu68806mTJnCM888w1VXXcVV1XdvgYkTJ7LFFluw+eabL/c5adIkhg3zrzdJkhpq+PBipNUOO8BHPlKsqyS1br8d7roL7ryzWKZPh6lTi+VPf1p2jCFDYKONYIstXr6su66PIkrSIOK/8CQ1hYhg0qRJTJo0aenk8ABz5szhnnvuWZrQqnxOnz6dadOmMW3aNK655prljjVs2DA22WQTtthiCzbbbLOlx9144415xStewSqrrLKSr06SJAHLJ7WqzZixfELrrruKRNesWcsSW5ddtvw+o0cvm2h+442LZdKk4nODDYpJ5SVJLcO/1SU1tdGjR7PLLruwyy67LLf+hRde4MEHH+T+++/ngQceWPr5wAMPMG/ePO6//37uv//+To85ceLEpQmtSZMmseGGGzJ9+nRe9apXsdFGG9HW1rYyLk2SJFWsuSbsuWexVGTCM8/A/fe/fJkyBebMgX//u1g6Gjq0eJyxTGgN2Wgj1nvxRWLNNYt166xTlJEkDRgmsCQNSOPHj2fXXXdl1113XW79kiVLePLJJ5cmsB5++GGmTp3KlClTmDp1KrNnz146cuv6669fbt+vfOUrRAQTJkxggw02YIMNNmD99ddf7nODDTZg4sSJPqIoSVJ/iygSTeusA69//fLbFi4sRmXdfz88+GDxfcqUZaO1FixY9v2qqxgK7Azwve8V+w8bBuutV4zUWn/94rPjstZaPqIoSU3Ef4FJailDhgxZmmh605vetNy2zOT5559fmsyqJLamTJnC3XffzYwZM2hvb1+a4Lrpppu6PMfEiRNZf/31WXfddZk4cSITJkxgwoQJS79PnDiRtdde29FckiT1h+HDl82F1dGSJcWcWpWE1pQpLHnoIZ6/7TbWmD2beOopWLQIHn20WLo7x/rrF3NtTZxYLBMmLP85cWIxemzIkP67VkkSYAJL0iASEayxxhqsscYayz2S2N7eziWXXMJ+++3Hiy++yOOPP87jjz/OE0888bLPJ598kvb2dp588kmefPLJHs+55pprLpfUqiS61l57bdZaay3WWmst1lxzTdZaay3n5pIkqR6GDCmSTuuuC699LQCL29v55yWXsP/++9MWUSS4Hn+8WJ54Ytn3yjJ9ejHKa8qUYunO0KHFKLGOia211nr5suaa4M0tSVohJrAkqTRkyBDWWWcd1llnHXbeeedOyyxZsoRnnnlmaZKrMpl89WdlWbx4MTNmzGDGjBnceeedPZ5/1VVXfVlSq6vfq6++OuPGjWOId3wlSarNsGHFyKr114fdd++8zMKFMG1akcyaNq1Ypk9/+fdnn4XFi+Gpp4qlN1ZbrfPkVnWSa621YPXVi2XsWB9llCRMYElSTYYMGbJ0FFXHieWrLVmyhOeee+5lya3K44nPPvssM2bMWPq5aNEiXnrpJV566SWmTp3aq1gignHjxrH66qszfvz4Hj+rv48ePZqwMyxJUueGD4eNNiqW7rS3F0mszpJczz67/PLcc8XjjS++WCwPPti7WIYOLZJelYTW+PG9+z5+PIwY0ceKkKTmYQJLkvrBkCFDlo6W2nbbbbstm5m8+OKLSxNa1cmtrn7PnTt36X4vvvhizfG1tbUxfvx4xo0bx7hx4xg7dmyvv48dO5ZRo0aRmStYO5IktYi2tmWPK/Zk8WJ44YWXJ7Y6W2bMKMrOm1fs99xzxVKrUaNg3LhiFNe4ccsvPa0bNYohCxcWb4OUpCZgAkuSGiwilo6Q2myzzXq1z8KFC3nhhRd4/vnna/5ctGgR7e3tPPPMMzzzzDMrHPeQIUMYM2bMcgmusWPHsuqqq7LqqqsyZsyYpd+7W1dZP2rUKEeFSZJa19ChxeOBa64JW23Vu33mzSsSWS+8AM8/Xyydfe+47oUXisTT3LnFMm1azeG2AW8Hsq1t+STXmDGw6qrLllp/r7KKj0RKWiEmsCRpABo+fPjS+bpqkZnMmTNnaUJr5syZzJo1a7nP3nxfvHgxS5YsWbq+HiKC0aNHd5vwGjVq1AovI0aMMEEmSRpYRo4slt6M8Kq2ZAnMnFkksmbNKr5Xlo6/u1iXs2cTmUR7ezEibMaM+lzT0KFdJ7hGjy5GjVU+O1u62zZqVPH4p6SWZAJLkgaRiFiaDNpwww1X6BiZycyZM/nDH/7AzjvvzLx585YmsmbPnr10Lq/Olq62ZyaZufR3f4iIXie7Ro4cyciRI1lllVWWLh1/92b9sGHDTJpJkla+IUOWzYO1ghYtWMDk3/+efXfbjba5c5cluV56CWbPLj4rS/XvrrbNnVscePHiZYmy/jBsWPcJrkoCbOTIYjRYZenud09lhw7tn2uRtBwTWJKkmlRGSq2++upsueWWtPXxdeCZybx583pMcs2ePZt58+Yxd+7cmpaFCxcuPc+cOXOYM2dOPaqhV4YMGbJcUmvx4sWsscYaPSbDRowYsXQZPnx4p5+9XVf96VsrJUm9NmQIi0aNgg02KOb66qvFi2HOnO4TX5VHHufMWfa949LZtjlzilFnAIsWFcm2WbP6HnNvtbXBKqswbJVV2CeTYePHL0tydZb8GjFi2TJ8eG3feyo3bJiPaKplmcCSJDVU9ciotddeu+7HX7RoUc2Jrzlz5jB//vyXLfPmzet0ffW2BQsWLD33kiVLlh6z4umnn677NfbWsGHDak58VX9va2tj+PDhtLW1Lfe942d323pbZtgwuyiS1FKGDi0miR87tv7HzizeCNlT8qt6+7x5MH/+sqX6d2+2LVq07Pzt7dDeTsyezSio3+OWKyJixZJjw4cXS1tb/T57U0aqgb1DSVJLGzZsGGPGjGHMmDEr5XxLlixh4cKFL0tszZ49m6uvvpodd9xxaVKtu2TYwoULWbBgQZefvdnW3t6+XGyLFi1i0aJFK3UUWl8MGzaMt73tbey///6NDkWS1MwiliVgVltt5Zxz0aLlk1zz59M+axbXX3UVr9lpJ4a1t3efBFu4EBYsWPbZl+/Vb4rMXHaeAWDY0KG8YcMNwbZevWACS5KkOqp+bLBae3s706ZNY8899+zzY5e9tWTJEtrb22tOfHW1rb29nfb29qXJsY6fK7pt4cKFLKq+k13qbJ0kSU1h2LBlk89XtLfz4uOPk6997cobXZRZJNP6kghbsKAYRbZw4fKfna3ry2c5rUO1WLyYqE7ASd0wgSVJUosaMmTI0scAm11msmjRouWSWnPnzuWf//xno0OTJKl5RSx7XG/06EZH073MYi60qsRW+5w53HDVVezd6Ng0IJjAkiRJDRcRS+fGqmhvb2fV6jvbkiRp4IooRq4NG1ZMZg8wbhwLVl+9sXFpwPB1RE0qIl4fEX+OiKciIiPiwF7ss2dE3BYRCyLioYg4rJMyR0XEIxExPyJujIhdO2xfJSJ+GBHPRcRLEfG7iFinflcmSZKq9dQ2d1L+4Ii4ryx/Z0Ts32F7RMQJETEtIuZFxJURsVmHMqtHxK8jYlZEvBgRZ0WE2UJJktS0TGA1r9HAHcBRvSkcEZOAvwJXA9sDpwA/i4g3V5V5L3Ay8A1gx/L4l0dE9Wu/vg+8HTgYeAOwLvD7vl2KJEnqTC/b5uryewDnA2cBOwAXAxdHxNZVxY4HPgUcCewGzCmPWT0x26+BVwH7AG8DXg+cUbcLkyRJqjMTWE0qMy/NzC9n5h96ucuRwNTMPDYz783M04HfAv9dVeYzwJmZeXZm3lPuMxc4AiAixgEfBj6TmVdl5q3A4cAeEfHqOl2aJElaptu2uROfBi7LzO+W7f1XgNuAo6EYfQUcA/xPZv4xM/8DHEpxQ+rAssxWwH7ARzLzxsy8Dvgv4H0RsW4/XackSVKfmMBqHbsDV3ZYd3m5nogYDuxUXSYzl5S/dy9X7QS0dShzH/BYVRlJklQHvWybO+q2vQcmARM6HHMmcGNVmd2BFzPzlqpjXAksoRixJUmS1HScxL11TACe7rDuaWBsRIwExgNDuyizZdUxFmbmi52UmdDViSNiBFD9iqsxwNLXpvdFZf++Hmcwsc5qZ53VzjqrnXVWu3rWWZPW+5r03DZ31FV7P6FqO70o80z1xsxcFBHP00V7b1vfXKyz2llntbPOamed1W4QtPWqIxNYqocvAF/ruHLy5MmMGjWqLie44oor6nKcwcQ6q511VjvrrHbWWe3qUWdz586tQySDmm19E7LOamed1c46q511VjvbevWGCazWMR3o+LbAdYBZmTkvIhYDi7soM73qGMMjYrUOo7Cqy3TmRIoJaCvGAE/su+++jB07trar6KC9vZ0rrriCffbZZ7lXq6tr1lntrLPaWWe1s85qV886mzVrVp2iqqsZ9Nw2d9RVe1/dllfWTetQ5vaqMstNEh8Rw4DVuzmvbX0Tsc5qZ53VzjqrnXVWu0HQ1quOTGC1jhuA/Tus26dcT2YujIhbgb0p3lhERAwpf59elr8VaC/X/a4sswWwYeU4ncnMBcCCyu9i/lhoa2ur21/c9TzWYGGd1c46q511VjvrrHb1qLNmrPNets0d3VBuP6Vq3dL2HphKkYTamzJhFRFjKea2+nHVMVaLiJ3KF7YA7EUxN+qNXcRqW9+ErLPaWWe1s85qZ53VrlXbetWXCawmFRGrAptWrZoUEdsDz2fmYxFxIrBeZh5abv8JcHREnAT8nKIj+h7grVXHOBk4NyJuAW6ieEvRaOBsKCZ5jYizgJPLeTBmAacBN2Tmv/rnSiVJGtS6bZsj4hfAk5n5hbL8qcDfI+JY4K/A+4CdgY8BZGZGxCnAlyPiQYqE1jeBpyiTZJl5b0RcBpwZEUdSvMDldOA3mflUf1+wJEnSijCB1bx2Bq6u+l0Ztn8ucBgwkWJkFACZOTUi3gp8n+IV209QvB778qoyF0TEWsAJFJO03g7sl5nVE73+N8VbiH5HMVnr5cAnV+QC6jGEs729nblz5zJr1iwz6r1kndXOOquddVY766x29ayzZn2soBdt84YU7XKl/PURcQjwP8C3gAeBAzPzrqrDnkSRBDsDWA24rjzm/KoyH6BIWv2NZe3+p2qN37a+Mayz2llntbPOamed1W4wtPWqn8jMRsegFhMR61Ek0CRJajbrZ+aTjQ5ioLOtlyQ1Mdv6FmUCS3UXxcQY6wKz63C4MRQd5PXrdLzBwDqrnXVWO+usdtZZ7epdZ2OAp9LOT5/Z1jecdVY766x21lntrLPa2dar13yEUHVX/mVRl4x3ZZJYYHZmOia0F6yz2llntbPOamed1a4f6sx6rxPb+sayzmpnndXOOquddVY723rVYkijA5AkSZIkSZK6YwJLkiRJkiRJTc0ElprdAuAb5ad6xzqrnXVWO+usdtZZ7ayzwcE/59pZZ7WzzmpnndXOOquddaZecxJ3SZIkSZIkNTVHYEmSJEmSJKmpmcCSJEmSJElSUzOBJUmSJEmSpKZmAktNLSKOiohHImJ+RNwYEbs2OqZmERGvj4g/R8RTEZERcWCH7RERJ0TEtIiYFxFXRsRmDQq34SLiCxFxc0TMjohnIuLiiNiiQ5lVIuKHEfFcRLwUEb+LiHUaFXOjRcQnIuI/ETGrXG6IiLdUbbe+ehARny///zylap31ViUivl7WUfVyX9V266vF2dZ3zba+Nrb1tbOt7zvb+p7Z1qteTGCpaUXEe4GTKd5KsSNwB3B5RKzd0MCax2iKOjmqi+3HA58CjgR2A+ZQ1N8qKye8pvMG4IfAq4F9gDZgckSMrirzfeDtwMFl+XWB36/kOJvJE8DngZ2AnYGrgD9GxKvK7dZXNyJiF+DjwH86bLLeXu5uYGLV8tqqbdZXC7Ot75FtfW1s62tnW98HtvU1sa1X32Wmi0tTLsCNwOlVv4cATwKfb3RszbYACRxY9TuAacBxVevGAfOB9zU63mZYgLXKent9Vf0sBA6qKrNlWebVjY63WRbgeeDD1leP9bQq8ADwJuAa4JRyvfX28rr6OnB7F9usrxZfbOtrqivb+trrzLZ+xerNtr539WRb3/u6sq13qcviCCw1pYgYTnEn6MrKusxcUv7evVFxDSCTgAksX38zKf6hYP0VxpWfz5efO1Hcqa2us/uAx7DOiIihEfE+itEAN2B99eSHwF8z88oO6623zm1WPiI1JSJ+HREbluutrxZmW99ntvU9s62vgW19zWzra2Nbrz4b1ugApC6sCQwFnu6w/mmKjLy6N6H87Kz+JjDIRcQQ4BTgn5l5V7l6ArAwM1/sUHxQ11lEbEPRiV0FeAl4Z2beExHbY311quz87wjs0slm/zt7uRuBw4D7KR4p+BpwbURsjfXV6mzr+8a2vhu29b1nW1872/qa2darLkxgSRqMfghszfLP3qtz9wPbU9zFPgg4NyLe0NCImlhEbACcCuyTmfMbHc9AkJmXVv38T0TcCDwKvAeY15ioJLUA2/res62vgW197WzrVS8+QqhmNQNYDHR8+8Q6wPSVH86AU6kj66+DiDgdeBvwxsx8omrTdGB4RKzWYZdBXWeZuTAzH8rMWzPzCxSTCX8a66srOwFrA7dFxKKIWEQxGemnyu9PY711q7wD+wCwKf531ups6/vGtr4LtvW1sa2vmW19H9nWa0WZwFJTysyFwK3A3pV15VDwvSmGOKt7Uyn+wq+uv7EUbygalPVXvmr8dOCdwF6ZObVDkVuBdpavsy2ADRmkddaFIcAIrK+u/A3YhuJOdmW5Bfh11XfrrRsRsSqwCcXk1P531sJs6/vMtr4D2/q6sa3vnm19H9nWa0X5CKGa2ckUQ5hvAW4CjqGYVPLsRgbVLMq/+DetWjWpnKvg+cx8LCJOAb4cEQ9SdHK/CTwFXLySQ20WPwQOAd4BzI6IyjP1MzNzXmbOjIizgJMj4nlgFnAacENm/qsxITdWRJwIXEoxieYYivrbE3iz9dW5zJwN3FW9LiLmAM9V5mCx3pYXEf8H/JniUYJ1gW9QjMo53//OBgXb+m7Y1tfMtr5GtvW1s62vnW296sUElppWZl4QEWsBJ1BM4Hc7sF9mdpysdLDaGbi66vfJ5ee5FJMknkTxj4AzgNWA6yjqb7A+q/+J8vOaDusPB84pv/83sAT4HcWdx8uBT66E2JrV2sAvKCbbnAn8h6JDe0W53fpaMdbb8tYHzgfWAJ6l+Lvq1Zn5bLnd+mphtvU9sq2vjW197Wzr+4f1tjzbetVFZGajY5AkSZIkSZK65BxYkiRJkiRJamomsCRJkiRJktTUTGBJkiRJkiSpqZnAkiRJkiRJUlMzgSVJkiRJkqSmZgJLkiRJkiRJTc0EliRJkiRJkpqaCSxJkiRJkiQ1NRNYkiRJkiRJamomsCRJkiRJktTUTGBJkiRJkiSpqZnAkkREHBYRGRGvaMC5vx4RubLP21sRsUtEXB8Rc8o62r7RMdWqlj/fVrje/hIRx0fEfRGxUtvOiDgyIh6LiBEr87ySWpttf9daoS207e+9RrXv5blt46UamMCSmkhVZ6Or5dWNjnEwiYg24CJgdeC/gQ8CjzY0qH60Mq43IvYo/+GyWj2P28M5Dy///9m1w/qIiAfLbTt22DY0Ip6IiKvK32OBzwHfycwlVeVGRMR3IuKpiJgXETdGxD51voRzgOHAx+t8XElNwLa/udj2t0bb31tdte/lNtt4qckMa3QAkjr1VWBqJ+sfWtmBDHKbABsBH83MnzU6mJVgZVzvHsDXKDpsL/bTOTqaWX6O7bB+X2DTLrYdAKwHHFP+PoKizTy/Q7lzgIOAU4AHgcOASyLijZl5Xd/CLmTm/Ig4F/hMRJyWmU07akFSn9j2Nwfb/vprRNvfW12172AbLzUdE1hSc7o0M29pdBBi7fLzxXodMCJGZ+aceh2vzup+vStDL+q0qwTWJ4B7gFd2se0p4OLy9+HAnzJzftV5dwXeB3w2M/+vXPcL4C7gJIoOe59FxCrAhcDxwBuBq+pxXElNx7a/Odj2DwB1rNOXte/l8fu9jY+IVcrz2sZLveQjhNIAExEHlY8UvKGTbR8vt21dtW6HiLg0ImZFxEsR8beeHkeo5RwRsV5E/Dwino6IBRFxd0Qc0cVxXxsRN0fE/Ih4OCJqGi5dnuuscij3goiYGhE/jojhtVxvOYw9I2LTiDgnIl6MiJkRcXZEjCrLnAP8vdzlorL8NSt4nldGxHkR8QJwXYdtm0fEr8rzPxsR34zCBhHxx/L40yPi2A7H3igifhQR95fD2p+LiItiBecy6ep6az1Pd39GEfF14Ltl0amx7PGYV9SjTrtRSWCNqTrOBsDbgJOBRR22bQq8CTgzMxdFxCRgW+DKDsc9CFgMnFFZUXZEzwJ2L89Rk4i4Iop5SF4XEX+PiHnAqZl5K/A88I5ajylp4KulXS7X9fj3aV/OEbb9tv3LH6fp2v6I+Ej539w/I2KjqvUREVdHxIyI2I3O23eocxsfXbTv5XFt46VecgSW1JzGRcSaHdZlZj4H/BV4CXgPyzodFe8F7s7MuwAi4lXAtcAsirtF7RTP2F8TEW/IzBu7OH+vzhER6wD/AhI4HXgWeAtwVkSMzcxTKjtFxDbA5LLM1yn+/vkG8HSPtVHsvy5wE7AaRWfiPopHvA4CRgELV+B6L6R4XOMLwI7AR4BnKOZC+CnwJPBF4AfAzZVYV+A8F1EMPf8iEB22XQDcC3weeCvwZYpOzMcp7sJ9DvgA8H8RcXNm/qPcbxeKu3+/AZ4AXkExauiaiHhlZs7toUo76up6e32eXvwZ/R7YHHg/xTwbM8pdn61znXY0q/ysHmX1sXL9eRQd6+ptR7J8p7Vyl/W2DsfdAXggM2d1WH9T+bk98HgPsXW0LcVd8IuBM8v4Ko8P3Qa8psbjSRo4bPs7sO237e8Qby1t/80U7fsXgeOA/yrXHwXsCRzCsmkEOrbvUP82vrv2vRKDbbzUk8x0cXFpkoXi2frsYplfVe48ig7G0Kp1Eyj+0f2VqnV/ABYAG1etm0jRUfh7J+d9RS3nAH5G8ZjVGh2u43yKRnpkh1jmARtWrduKYvRL9qJuzi3PvXMn26LG6/16eb1ndTjO74EZVb/3LMsd1KFcrec5r5OYK9t+WrVuKEVnaAnwuar1qwFzgXOq1o3s5JivLo/5wS7+u3pFx306lHvZ9dZ4nt78GR3XWSz1qNNurmudcp8vlb/bgGnA/5W/H63UN7AK8Bzw26r9v1nuv2qH494F/K2T872yLP/x3sZY7rd2ud9sYMtOtv8UmFvLMV1cXJp/wba/u7qx7bft77ZOe/Hf0OXA9eX3jSmStH8of3favpfb6tbG00P7XpaxjXdx6cXiI4RSczoK2KfD8paq7RdQNIZ7Vq07iOKx4AugeIsaxSTVF2fmlEqhzJxG0UF9bRRvXulKt+eIiADeDfy5OF2sWVkoOgvjKO5sVmJ5cxnLY1Wx3FuW7VYUrzU+EPhzdjI/SGbmCl7vTzr8vhZYo7t6qdN5qi2dMDUzFwO3UNxVPKtq/YvA/RQdr8q6eVUxtUXEGhR38l6krPd66O15evNn1NU5+qFOO+o4B9a7KJJaPy5/z67a9h6KNzH9sGr/NYBFmflSh+OOpOh4dzS/anstti0/v5WZ93Wy/QVgZJSPukhqObb9VWz7bfv72PZX3A5sXcb6c4p2+xPltq7ad6hvG99T+w628VKvmMCSmtNNmXllh+Xqqu2XUfyj/L1V694L3J6ZD5S/16IYun1/J8e/l+L//+6e3+/pHGtR3B38GMWjAdXL2WX5ysSga1E09A92cp7O4utoLYoEw109lKn1eh/r8PuF8nN8nc8ztZvjdYxhJsUd9xmdrF8aV0SMjIgTIuJxig7WDIq6X43iHxB1UcN5evNn1JV61+lyspizYiHL5rn6BHBZZj5c/p7dYdu9Hf5/68o8YEQn61ep2l6LbcrPC7rYXnlcost/EEga0Gz7l2fbb9u/wm1/lbso2vjvAm8APp2Z03uxXz3b+J7ad7CNl3rFObCkASgzF0TExcA7I+KTFKNJXkPxnP/KOkclAf4riuHjnflPveLpJ4u7WN/TvAq16q6T01kMvYnrNIo355wC3EDRyU2K+SrqeXNiZZ2nVrUmh2YBYyPilRQd2LdWbZtdbtuO4hGJT3XY9zlgWESMyczZVeunUczz0dHE8vOpGmPcFphWfSe6g/EUjxfUeu2SWoBtf93Y9vesVdp+WJZc+wzwl8z8VdW2rtp3qG8b31P7DrbxUq+YwJIGrguADwF7U8wnESx/Z+dZirkTtuhk3y0p5lroafLJ7s7xLMU//IdmZmdvb6n2LEWnY7NOtnUWX2f7zwK27qFMX6+3N1bWeXpyEHBuZi59Q1FErEJxd7QR5+nNnxF0fmdxZdTpTIq7xJ8AplCMMqiYXbVtDi//R1lluP8klv+H2e3AG8tJi6sned2tansttgXu6Gb7JIq70pIGL9v+l5ex7bft705lhNeLFBPEV+uqfYf6tvE9te+VGGzjpR74CKE0cF1J8caa95bLTZm5dGh1OafCZOAdUfXa4/LtQYcA1+XL36zS63OUx/8d8O6oenV31XnW6hDL5cCBEbFhVZmtKObH6FZmLqF4a8vbI2LnTs4VdbreHq2s8/TCYl5+t/i/KCaDXenn6c2fUfl1Tvm5WtW+K6NOZ1LcNT0U+HEZb8VsirusHwB+3cm5big/O17Xbynq4WNVMY+guGt9Y2b2uuNdzgXySrrv4O4IXN/bY0pqSbb9y45j21+w7e/eR8vPP2Vmx1FTXbXvUKc2vpftO9jGS73iCCypOb0lIrbsZP31leHHmdkeEb8H3geMpnjDS0dfppgE9rqI+BHFW38+TvFM//E9BdGLc3weeCNwY0ScCdxDMQH2jsCbyu8VXwP2A64tYxlG0Rm6m2WTW3bnixSTff49Is6guEs1ETgYeC3FnbU+XW8NVtZ5uvMX4IMRMZOi3nenqPPnGnie3vwZ3VqW/d+I+A3FK7P/TP/X6UyK/1bnUUziWm02sGv5/Ucdd8zMKRFxF8V1/7xq/Y0RcRFwYkSsTTHB7YcoXjf+4Y7HiYikeKvSnp3EtxnFvBqddnAjYieK/5/+2OUVShrobPtfzrZ/ebb9NYiITYBvlT9f9v9WV+17ua3XbXxf2vdyf9t4qZdMYEnN6YQu1h9O8fhTxQXARyiGZl/YsXBm3h0RrwNOBL5AMeryRuD/ZeaNvYyly3Nk5tMRsSvwVYo3u32SonNzN/C5DmX/ExFvBk4ur+8Jio7tRHrRic3MJyNiN4pXHn+A4pGvJ4FLKYag1+t6e7SyztODT1PcIf0ARcfonxQdsB7f7NRf5+nln9HNEfEV4EiKf9QMASathDqtvInw/Mx8vsO2yrwX12dmVx3MnwMnRMTIDvNTHEpxvR+kmL/iP8DbMvMf1TtHxKrl12ldHL8ywWtXc8ccTDHp71VdbJc08Nn2v/xabPuXZ9vfS+Xor7MoJqG/ADi4HLXX8XHGrtp36EUbX4f2HWzjpV6Lbt5uKkmSgIgYR/EPyOMz86yeyney//4Ud7S3y8w7a9x3BPAI8O3MPLXWc0uSNNhExFHA6RRJqHbgfGCTjhOpN7J9L/e3jZdq4BxYkiT1IDNnAicBn42IFWk73wj8ZkU6txSjL9qBn6zAvpIkDSrlnFrfBv6cmb8EKm3vjh3LNrh9B9t4qSaOwJIkSZIkDXjlo4NXAjsAr8rMaRExDHgBeAr4HsULW+Z0cxhJTcoRWJIkSZKkVvAxYC/g05k5DSAzF1G8PGAU8ANgYePCk9QXjsCSJEmSJElSU3MEliRJkiRJkpqaCSxJkiRJkiQ1tWGNDkCtp5w8cV1gdqNjkSSpyhjgqXT+hD6zrZckNSnb+hZmAkv9YV3giUYHIUlSJ9YHnmx0EC3Atl6S1Kxs61uUCawmFRGvBz4L7ARMBN6ZmRf3sM+ewMnAq4DHgf/JzHM6lDmqPO4E4A7gvzLzpqrtq1C8XvZ9wAjgcuCTmfl0DeHPBnj88ccZO3ZsDbu9XHt7O5MnT2bfffelra2tT8caLKyz2llntbPOamed1a6edTZr1iw22GADGGAjhhrVH+gF2/oGss5qZ53VzjqrnXVWO9t61cIEVvMaTdGh/Dnw+54KR8Qk4K/AT4APAHsDP4uIaZl5eVnmvRQd2iOBG4FjgMsjYovMfKY81PeBtwIHAzOB08vzv6bWCxg7dmxdOrWjRo1i7NixNgK9ZJ3VzjqrnXVWO+usdtYZ0Lj+QK/Y1jeGdVY766x21lntrLPaWWeqhQmsJpWZlwKXAhTTTPToSGBqZh5b/r43Il4L/DfFKCqAzwBnZubZ5XGPpEhWHQF8OyLGAR8GDsnMq8oyh5fHenVm/qsuFyfp/7d35/Fy1vXd/1+fs+dkXwh7wk5AENlFRW1ZVKxLtcX9dvlVxVqV4l0qN+5a9dYWseJS6S2itQqVFqtEIOBCFQRkCwIhhIQEwhKSkJwkZ5tzzvf3xzVDJodzkkzOTK45c17Px+N6XDPXXDPzOZ9Sryvv+V7fS5J2Sh7nA9WrXpIkqXoMsBrHKcANw7ZdB1wMEBFtZJcffLH0YkppKCJuKL6X4uut5Z+TUloSEauK+4wYYEVEO9nlhiVTIUvTC4XCrv9Fxc8oX2vH7Fnl7Fnl7Fnl7FnlqtmzCdT3apwPPIfH+vpizypnzypnzypnzyrnsV6VMMBqHHsBw+epegqYFhGTgJlA8yj7LCj7jP6U0oYR9tlrO999AfCp4Ruvv/56Ojs7d6r4HVm0aFFVPmcisWeVs2eVs2eVs2eVq0bPuru7q1DJuFCN84GReKyvQ/ascvascvascvasch7rtTMMsFQNXySbS6NkKvDYmWeeWZV5MRYtWsQZZ5zhNdE7yZ5Vzp5Vzp5Vzp5Vrpo96+rqqlJVE5bH+jpizypnzypnzypnzyrnsV6VMMBqHE8Cew7btifQlVLqiYhBYHCUfZ4s+4y2iJgxbBRW+T7PkVLqA/pKz0tzdLS2tlbtf7ir+VkThT2rnD2rnD2rnD2rXDV6NoF6Xo3zgefwWF+f7Fnl7Fnl7Fnl7FnlPNZrZzTlXYCq5hayOw2VO6O4nZRSP3BH+T4R0VR8fktx0x1AYdg+hwPzyvaRJEn1qxrnA5IkSXXHEVh1KiKmAIeUbTowIl4ArE8prYqILwL7ppT+V/H1bwN/ExFfJrvV9p8CZ5PdVajkIuDyiPgDcBvZbbMnA5cBpJQ2RsT/Ay6KiPVAF/B14BbvQChJ0u6Xx/mAJElSPTLAql8nAL8qe16ad+Jy4F3A3mQjowBIKa2IiFcDXwU+AjwG/FVK6bqyfa6IiD2Az5JN8no38MqUUvlErn8LDAFXkd1t6Drgr6v5h0mSpJ2W1/mAJElSXTHAqlMppV8DsZ3X3zXKe47dwedeAlyyndd7gQ8WF0mSdouUoFDIlv7+bL1lC2zePLHns8jrfECSpFoYTIn+oSH6i+vu/n7Wx6iHOWkbBliSJDWolGBgIAuE+vpGXkZ7rRQi9fdv+3j4utLXRtu/UBjpL2jl9a8/lLPP3t2dkyRp/Egp0Z8SvUND9JUvKT3neSk8KpSFSIXi+599PMbXRtyv+DyNUP8Bkyfz9t3eNY1HBliSJFXZ0FAWAvX2Qk9Ptu7qgmXLpvO73wUDA1u3l9a9vTsOlXYUOI20PY10pjhOtLYmUvJXWUlSfRoqhkY9Q0P0DA6yqb+fR5qauH3TJgaamugZGqK3uOwoWCp/3rsT+5Q/7x/HB/uWiNGHGUvDGGBJkhre4GAWFHV3Z5eldXeP/LinZ9tAaXjItLPrvr6RqmgFXr57//BhIqC9fdulrW3kbaXtra1bn5cej2Xbzr6npQUGBgZYuPA+YH6ufZMkjQ+DKdEzOMiWoSG2DA7SXXzcPThIdzFkejZwKgZFpfBpm+elxyNsL/+MEYOjqVNh8eLd/8eXaY2gvamJjqYm2ouPS0tbBG3FdWv54yq81hqx7X47eK01gsGBARYuXJhrvzR+GGBJknKXUjaCaNMm2Lz5uesdBU87etzbm9/f1tQEkybBpEmJlHqZObODSZOCSZOgo4Nn1x0dI4dJo4VM2wufRtu3pSULsSRJykvf0BCbBwfZNDDApsFBtpSFTMMfdxeDqPLHI20rPe4dGsrt72qJYFJTE02FAtMnTWJSUxOTmptHDJFG2raj5zv7nramJprG0cF+MO8CNK4YYEmSdsnAQLB2bRYQbdw4cvBUvt7RawMDu6fuzs5tl8mTtz7OgiaeEy5Vui5/3Fqcg7xQGGDhwus566yzaG2d2BOTS5LGj4GU2FQosKksdCotm0fYVnq+eYRtmwYHKeymy906m5qY3Ny8dV0MkyYVl2cfNzdv+7z0eITtk0b5jI6mJlqamigUCixcuJCzTj3VY71UAwZYkjQBFQpZ6FRaNmzY9vmOlxZ6el5bk9omTYIpU7IR+FOmZEspZCoPmyp5XFp3dGQjoiRJmghSSmweHGTjwAAbBgbYWP54+Hpw8Dnb102bRt/NN9ekto6mJqY2NzO5uZnJZSHT8MedZfuM9LgUUJU/ntTURIyjUUiSdo4BliSNUyllcy6tX1/5smXLWL9960nhlCkwbVoWOJVCp51Zj7Rt8uTsMjdJkrTVwNAQzwwMsK5QYP3AAOsLBdYV1+tL28seP1MMoboGBsZ2iVZZCNQewdSWFqY2NzOluZmppaW4rbRMGWFb+fPSe1v8RUlShfxngiTViUIB1q6Fp5/OljVrtj4uX8qDqJEnC995U6bA9OmVL52dBW67bRFvfOMZTJrkEHlJknbWUEpsGBhgTX8/TxcKrCkUWNPfz5pCgaf7+1k3QiDVNTi2mYJaIpjR0sL05uZsXVxmlK+HvTYFuPO3v+X1p53GrI4O2gycJOXMAEuSaqi7G554YttlzZqRw6lnntm172hthVmzKltmzMhGTe3qaKdCAR54oOBoKUmSgMLQEE/19/N4fz9P9vdvDaTKwqlSYPV0ocDALs4DNaOlhVktLcxqbWV2cT2rpYXZra3bbJs5LJzalUvqCoUCTwwNMbu1lVbDK0l1wH96SNIu2LQJVq9+bjhVvjz+OHR1Vfa5TU0wezbMnQt77JEt5Y/32CN7vTyMmjzZO8tJklQLhaEhniwGU0/09WXr/n4e7+vbZv10oUClkdT05mbmtrWxR2src9vamNvayh6trcwphVHFcKoUTM1oafGyO0kTmgGWJA3T1wePPQaPPgqrVmXr4cuGDTv/eR0dsPfeW5c999waSg0PqmbOhObmmv1pkiSpaCglnuzv59G+Plb19rKqr49Hi+vS86cLhZ3+vJYI9mprY++2NvYsBVPFcKo8pJrb1sac1lbaDaMkqSIGWJImnP5+WLkSli/PlmXLmrj55hP57Gebeeyx7PK+nTFt2rbB1GjL9OmOkJIkaXcbTInVfX0s7+lh6ZYtLGpv58qlS1nd38+qvj4e6+ujsBOX8rUWg6l92trYu709W7e1sU97+zbrOa2tNHnAl6SaMcCS1JA2bIClS+Hhh7cGVaXl0UezO/ht1Qzss837Ozpg3jzYf/+ty/DnU6fuxj9IkiQ9R/fgICt6e3m4p4eHe3pYXvb4kd5e+ssP+B0d2aSTZZqAfdvbmdfezv4dHcxrb2decb1/ezv7trcz22BKkuqCAZakcWtgAB55BB58MFuWLNn6+Kmntv/ezk446KBsOeCAQbq77+dVrzqCAw9sYf/9s3mmPFeVJCl/qXip35Lu7ucsq3ZwO96WCA7o6OCg9naannqKUw87jAM7O9m/GFTt09bmvFKSNE4YYEmqe4OD2cipxYvh3nuzZckSeOih7G54o9l7bzj00K1BVfkyd+7WgKpQGGLhwuWcddYCWlt3z98kSZK2lVLisb4+Fm/ZwuLNm3mgLKjaNDg46vumNzdz8KRJHDxpEgd1dGTrSZM4uKOD/drbaWlqolAosHDlSs7abz9aPdhL0rhkgCWprqxdC/fcszWoWrwY7rsPenpG3r+jIwupFiyAww/PlgUL4LDDsjmqJElS/dk8MMB93d0s3rz52cBq8ZYtbBgYGHH/JuCgSZNY0Nm5zXL4pEnMbm0lHDYtSQ3PAEtSbp55Bu64A/7wh63LypUj79vRAc97Hhx9dLYceWQWVM2bB478lySpfm0ZHOSuTZu4vbj8YdMmlvX0MNL06S0RLOjs5OjJkzlq8uRng6qDJ03yrn2SNMEZYEnaLfr6srDq5pvh9tuzsGr58pH3PfhgeP7zt4ZVRx8NhxwCzc27t2ZJklSZwtAQ92zezG1lYdX9W7YwNMK+e7e18fzJk3n+lCnPrhd0dtJmUCVJGoEBlqSaWLcuC6t+97tsuf32LMQa7uCD4YQTti7HHgvTp+/+eiVJUuU2DQzw+64ufrtxI7/duJHfd3XRPfTcuGrftjZOnDaNE6ZO5cSpUzl2yhT2aGvLoWJJ0nhlgCWpKp55Bn71K7jxxmz9wAPP3WePPeDFL4aTT4YTT4TjjoOZM3d/rZIkaddsHBjg1xs28MtnnuG3Gzdy9+bNzxldNbOlhZOnTePEYlh1wtSp7N3enku9kqTGYYAlaZd0d2cjq264IQut7rwT0rDJLBYsyAKr0nLooVvv/CdJkupf/9AQv+/q4oZnnmHRM89we1cXw+8HeEBHBy+ZPv3Z5YjOTpo84EuSqswAS9JOW7kSrrkGfv5z+OUvn3tJ4BFHwOmnw5/+KbzkJTBnTj51SpKkXfdYby8/W7eOn69bx683bHjOJYGHTprE6TNn8rIZM3jxtGns19GRU6WSpInEAEvSqIaG4NZbs8DqZz+De+/d9vX99oPTTtsaWu2zTz51SpKkXTeUEndu2sTP1q3jZ+vWcdfmzdu8vkdrK6fPnMnpM2dy2syZzDewkiTlwABL0jZSgt//Hq68Ev7jP2D16q2vNTXBi14Ef/Zn2XLkkV4SKEnSeJRS4vddXVyxZg0/efppVvf3P/taAKdMm8ZrZs/mVbNnc/TkyV4SKEnKnQGWJFKCO+6AH/84C61Wrdr62tSpcNZZWWD1qlfB7Nn51SlJknZdSok7Nm3iiqef5so1a1hVNhfAlOZmzpw5k9fOmcNZs2Z5h0BJUt0xwJImsKeegn/7N7jsMrjvvq3bp0yB170Ozj4bzjwTvFJAkqTx64m+Pr7/1FN894knWNrT8+z2Kc3NvG72bN40dy5nzppFe1NTjlVKkrR9BljSBDMwkM1p9d3vwsKFMFi8lVBHRxZavfnN8IpXwKRJ+dYpSZJ2Xf/QENesW8d3n3ySX6xb9+ydAyc1NfGaYmj1qlmzmNTcnGudkiTtLAMsaYJYuxYuvRS+9S149NGt208+Gd79bnjTm2DGjNzKkyRJVfBkXx/ffvxxvv344zxVKDy7/cXTpvGevffmL/fYg6kt/hNAkjT+ePSSGtzdd8PXvgY/+hGUprqYMycLrd79bjjiiFzLkyRJVXB7Vxdfe+wxrnz6aQopAbBXWxvv3HNP3r333hze2ZlzhZIkjY0BltSgfvvb4Mtfhmuv3brt+OPhQx/KRls5r5UkSeNbSolfrFvHP6xcye+6up7d/qJp0/jIfvvx53Pm0Oq8VpKkBmGAJTWQlOC664ILLngJDzyQ/b93UxP85V/CRz4CL3wheBdsSZLGt6GUuKWlhc/ccw93bdkCQFsEb547lw/tuy8nTJuWc4WSJFWfAZbUIG66CT72MbjllhZgNm1tiXe9Kzj/fDj44LyrkyRJY5VS4qdr13LhihXcP3kybNnC5KYmztlnHz66//7s3d6ed4mSJNWMAZY0zi1eDBdckN1REGDSpMQZZzzM1742nwMOaM23OEmSVBW/3bCB85cv55bipYKdKfGR/ffnvHnzmNPWlnN1kiTVngGWNE6tWZONuPre97JLB5ub4b3vhQsuGOCuu+5j333n512iJEkao6Xd3fzvhx/mZ+vWATCpqYmP7LMPR913H2fPn09rqz9WSZImBmd1lMaZgQH4+tfhsMPgssuy8OpNb4IHHoBvfQv23jvvCiVJ0lhtGRzkwuXLOer22/nZunU0A+fssw8Pn3wyn50/nyl5FyhJ0m7mCCxpHLn1Vnjf+7LLBgGOOw6+8Y1scnZJkjT+pZS4eu1azl22jFV9fQC8atYsvnrIIRze2QlAoVDIs0RJknJhgCWNA7298MlPwj/9EwwNwcyZ8IUvZJcMNjfnXZ0kSaqGNf39/PXSpVy1di0A89vb+dqhh/La2bMJbyMsSZrgDLCkOvf738O73w1LlmTP3/52+OpXYc6cfOuSJEnV85M1a/jAQw+xtlCgJYLz99+fC+fPp9NfqiRJAgywpLo1OAhf/CJ86lPZqKu99oLvfAde85q8K5MkSdXSNTDAOUuX8qM1awA4evJkLl+wgGOnTs25MkmS6osBllSHnnwyG2l1443Z87e9Df75n2HWrHzrkiRJ1XPXpk2cff/9LOvpoRm4YP58PjF/Pm1N3mdJkqThDLCkOvOrX8Gb3wxr1kBnJ3zzm/DOd+ZdlSRJqpaUEv/y+OOcu2wZfSkxr72dHx95JKdMn553aZIk1S0DLKmOfOtb8KEPZZcPHnUUXHklHHFE3lVJkqRq6Rsa4gNLl3LZk08C8JrZs/neggXMam3NuTJJkuqbAZZUBwoFOPfcbLQVwFvfCpdemo3AkiRJjeHp/n7ecN99/HbjRpqA/3vQQXx0//29w6AkSTvBAEvK2ebN8IY3wKJF2fMvfAE+9jHwXFaSpMZx/5YtvPree3mkt5dpzc1cceSRvHL27LzLkiRp3HCGyDoXER+MiEciojcibo2Ik7azb2tEfDIiHi7uf09EvHLYPlMj4uKIWBkRPRFxc0ScOGyf70VEGrZcW6u/cSJbuxZOOy0LryZPhquvhgsuMLySJG0rj/MBVc+tXV2cetddPNLby8EdHfz+uOMMryRJqpABVh2LiDcBFwGfAY4D7gGui4i5o7zl88D7gQ8BRwLfBv4rIo4t2+dfgTOAdwBHA9cDN0TEvsM+61pg77LlLdX4m7TVo4/CqafCbbdldxf85S/hda/LuypJUr3J+XxAY3TD+vWcdvfdrB8Y4OSpU7n1+OM5YvLkvMuSJGncMcCqb+cBl6aULksp3Q+cA3QD7xll/3cAX0gpLUwpLU8pfQtYCHwUICImAW8Ezk8p3ZRSWpZS+jSwDPjAsM/qSyk9WbY8U/0/b+JatQpe+lJYsgT22w9++1s4adTf0iVJE1ye5wMag5+uXcur772XLUNDnD5zJjcccwyznaxdkqRd4hxYdSoi2oDjgS+WtqWUhiLiBuCUUd7WDvQO29YDvKT4uAVo3sE+JS+PiDXAM8AvgY+nlNaNUmt78btLpgIUCgUKhcIope6c0vvH+jn1ZPVqOP30Fh55JDjkkMS11w4wb142kXs1NGLPas2eVc6eVc6eVa6aPRuvfa+D8wHtooXr1vGX991HISXeOGcOPzzySNqb/O1YkqRdZYBVv+aQnVw+NWz7U8CCUd5zHXBeRNwEPAycBryh+DmklDZFxC3AJyLigeJnvYXsBHhZ2edcC/wnsAI4GPgC8IuIOCWlNDjC914AfGr4xuuvv57OKt1Gb1FphvNx7pln2vn4x1/M6tVT2XPPLXzsY7/lj3/s5Y9/rP53NUrPdid7Vjl7Vjl7Vrlq9Ky7u7sKleQiz/OBZ/ljVWV+uWEDb7j/fgop8Zdz5nD5YYfRNDhIYXCk06jKNWLPas2eVc6eVc6eVc4fq1SJSCnlXYNGEBH7AKuBF6WUbinb/mXgZSmlk0d4zx7ApcBrgER20noD8J6U0qTiPgcD3wVeCgwCdwJLgeNTSkeMUstBxc86PaV04wivj3RS+9jatWuZNm1apX/6NgqFAosWLeKMM86gdZwPue/qgj/5kxbuvTeYNy9xww0DHHBA9b+nkXq2u9izytmzytmzylWzZ11dXcyZMwdgekqpqyoF7gb1cj4QEZ9mhB+r/v3f/71qP1Y1igeam/nU5Mn0R3ByocDfdXf7i7Ek7Qbd3d289a1vhXF2rNfO83hav9aSnVDuOWz7nsCTI70hpfQ08PqI6ABmA48DXwKWl+3zMPCyiJgMTEspPRERV5TvM8LnLo+ItcAhwHMCrJRSH9BXeh7FW+i1trZW7R9p1fysPBQK8Ja3wL33wp57wi9/GRx8cG3/nvHeszzYs8rZs8rZs8pVo2fjuOf1cj7wRbKJ5EumAo+deeaZ/lhVZmlPD+9ZvJj+gQFeOXMm/7FgQU0uG2yknu0u9qxy9qxy9qxy1f6xSo3NAKtOpZT6I+IOsmH/VwNERFPx+SU7eG8vsDoiWskmab1yhH22AFsiYibwCuD80T4vIvYjOwF+Ypf+mAkuJXj/+2HRIpg8Ga65Bg4+OO+qJEnjQb2cD/hj1Y493d/Pa++//9m7DV511FF0NjfX9DvHe8/yYM8qZ88qZ88qN8F/rNJOMsCqbxcBl0fEH4DbgHOBycBlABHxfWB1SumC4vOTgX2Bu4vrT5PdafLLpQ+MiFcAATxINqLqK8CSss+cQnaJwFVkv+weXHz/MrI5NVShL30JLrsMmprgyivh+OPzrkiSNM7s9vMBVaZncJDX/vGPLO/t5cCODv776KNrHl5JkjTRGGDVsZTSFcV5LD4L7EV2IvrKlFJpItd5wFDZWzqAzwMHAZvJbpn9jpTShrJ9ppNdBrAfsJ4sqLowpVSa8W4QeD7wTmAG2WUH1wOfKP76qgpcfz1ceGH2+BvfgLPOyrceSdL4k9P5gHZSSokPPvQQv+/qYmZLCwuPPpq5bW15lyVJUsMxwKpzKaVLGOUSgZTSy4c9/w1w5A4+70pGuISg7PUesksINEaPPJLNe5US/NVfwTnn5F2RJGm82t3nA9p5//rEE1z25JM0AVceeSQLJk/OuyRJkhpS9WeVlERfH7zxjbB+PZxwAnz963lXJEmSqu32ri7+5qGHAPiHAw/k9Fmzcq5IkqTGZYAl1cCFF8Kdd8Ls2XDVVdDRkXdFkiSpmjYODHD2/ffTnxKvnzOHv583L++SJElqaAZYUpXdeCP80z9ljy+7DDyflSSp8Xz4oYd4pDhp+/cWLHj2zoySJKk2DLCkKlq/Ht75zuzx+98Pr3lNvvVIkqTq+8maNXz/qadoAn5wxBFMb3FaWUmSas0AS6qiD38YVq+Gww7bOgpLkiQ1jif6+nj/0qUAfGzePF48fXrOFUmSNDEYYElVct118MMfQlMT/OAH4E2IJElqPB966CHWDwxw3JQpfOqAA/IuR5KkCcMAS6qC7m74wAeyxx/6EJx0Ur71SJKk6vvZ2rVctXYtzcB3FyygrclTaUmSdhePulIVfPazsGIF7L8/fO5zeVcjSZKqbfPAAH/z0EMAnLf//hwzZUrOFUmSNLEYYEljtGTJ1vmuvvENmDo133okSVL1ffqRR1jV18cBHR1eOihJUg4MsKQx+ru/g4GB7I6D3nVQkqTG81B3N19bvRqAbxx6KJObm3OuSJKkiccASxqDG26An/8cWlrgH/8x72okSVIt/P3y5QykxFmzZnHW7Nl5lyNJ0oRkgCXtosFBOO+87PEHPwiHHZZvPZIkqfp+s2ED/1WcuP0rBx+cdzmSJE1YBljSLvrBD+Dee2HmTPjkJ/OuRpIkVdtQSnx02TIA3rfPPhw5eXLOFUmSNHEZYEm7oFDYerfBCy6AWbPyrUeSJFXf1WvXcsfmzUxtbubTTtwuSVKuDLCkXfCDH8Dy5TB3Lvz1X+ddjSRJqrahlPjMI48AcO5++zG3rS3fgiRJmuAMsKQKFQrw+c9nj88/H7yaQJKkxvNfa9eyeMsWpjU3c+5+++VdjiRJE54BllSh738fVqyAPfeED3wg72okSVK1lY+++sh++zGrtTXfgiRJkgGWVImhIfjKV7LHf/d30NmZbz2SJKn6/nvtWu519JUkSXXFAEuqwC9+AQ8+CNOmwXvfm3c1kiSpFi567DEAPrDPPo6+kiSpThhgSRX4p3/K1u97XxZiSZKkxnJ7Vxf/s3EjLRF8yNFXkiTVDQMsaSfddRf86lfQ3Awf/nDe1UiSpFr4anH01ZvnzmXf9vacq5EkSSUGWNJOuvjibH322bD//rmWIkmSauDR3l6uXLMGgPMcfSVJUl0xwJJ2wvr1cMUV2eOPfCTfWiRJUm18+/HHGQRePmMGx06dmnc5kiSpjAGWtBP+7d+grw+OOQZOOinvaiRJUrUVhoa47MknAfjgPvvkXI0kSRrOAEvagZTgO9/JHr/vfRCRbz2SJKn6rlm3jif6+5nb2spr58zJuxxJkjSMAZa0A7fcAvfdB5Mmwdvelnc1kiSpFi594gkA3rXXXrQ1eYosSVK98egs7UBp9NWb3wzTp+dbiyRJqr5Vvb1cu349AH+19945VyNJkkZigCVtR3c3XHVV9vg978m3FkmSVBuXP/kkQ8CfzJjBoZ2deZcjSZJGYIAlbcfPfw6bN8MBB8CLX5x3NZIkqdpSSvzwqacAeOdee+VcjSRJGo0BlrQdP/xhtn7rW528XZKkRnT35s082NNDR1MTf+7k7ZIk1S0DLGkU69fDL36RPX7rW/OtRZIk1ca/r1kDwGtmz2ZaS0vO1UiSpNEYYEmjuOoqKBTg+c+H5z0v72okSVK1DaXEj4qXD7517tycq5EkSdtjgCWN4sc/ztaOvpIkqTH9z8aNrO7vZ3pzM6+aPTvvciRJ0nYYYEkjWL8efvOb7PFf/EW+tUiSpNr4j+Llg2/YYw/amzwtliSpnnmklkZwzTUwOAhHHw0HH5x3NZIkqdpSSvx03ToA3uDk7ZIk1T0DLGkEV1+drV//+jyrkCRJtXLn5s081tdHZ1MTp82cmXc5kiRpBwywpGF6euDaa7PHBliSJDWmn65dC8ArZ81iUnNzztVIkqQdMcCShrnxRujuhv33h2OPzbsaSZJUC6UA63VePihJ0rhggCUN89OfZuvXvQ4i8q1FkiRV34qeHhZv2UIz8GrvPihJ0rhggCWVSQmuuy57/OpX51uLJEmqjWuKk7efOmMGs1tbc65GkiTtDAMsqcyDD8Kjj0J7O7z0pXlXI0mSauH6Z54B4FWzZuVciSRJ2lkGWFKZ66/P1qeeCp2d+dYiSZKqr39oiF9t2ADAmd59UJKkccMASypTunzwzDPzrUOSJNXGLV1dbB4cZG5rK8+fMiXvciRJ0k4ywJKK+vrg17/OHr/iFbmWIkmSauT69esBOGPmTJq8W4skSeOGAVadi4gPRsQjEdEbEbdGxEnb2bc1Ij4ZEQ8X978nIl45bJ+pEXFxRKyMiJ6IuDkiThy2T0TEZyPiieI+N0TEobX6G+vFzTdDdzfstRccfXTe1UiSpFoozX91pvNfSZI0rhhg1bGIeBNwEfAZ4DjgHuC6iJg7yls+D7wf+BBwJPBt4L8i4tiyff4VOAN4B3A0cD1wQ0TsW7bP+cCHgXOAk4Etxe/tqNKfVpduvDFbn3EG+IOsJKme1OAHreaI+FxErCj+WPVwRHwiorGPgGv7+7lj0yYgG4ElSZLGDwOs+nYecGlK6bKU0v1kgVI38J5R9n8H8IWU0sKU0vKU0reAhcBHASJiEvBG4PyU0k0ppWUppU8Dy4APFPcJ4Fzg8ymln6aUFgP/C9gHeH1t/sz6cNNN2fplL8u3DkmSytXoB62/Jzv2/w1wRPH5+cX3NKybNm4kAc/r7GTv9va8y5EkSRVoybsAjSwi2oDjgS+WtqWUhiLiBuCUUd7WDvQO29YDvKT4uAVo3sE+BwJ7ATeUfe/GiLi1+L0/HqHW9uJ3l0wFKBQKFAqFUUrdOaX3j/VzdqS3F269tQUITjmlQI2/rqZ2V88aiT2rnD2rnD2rXDV7Ns77/uwPWgARcQ7warIftL40wv7vAP4hpbSw+PxbEXE62Q9aby9uexHw05TSNcXnj0TEW4BRR3Y1gv/ZuBGAl86YkW8hkiSpYgZY9WsOWdj01LDtTwELRnnPdcB5EXET8DBwGvCG4ueQUtoUEbcAn4iIB4qf9RayYGpZ8TP2Kvue4d+7FyO7APjU8I3XX389nZ2do7ylMosWLarK54zmvvtm0d9/KjNm9LJ06XU89FBNv263qHXPGpE9q5w9q5w9q1w1etbd3V2FSna/Gv2gBXAz8L6IOCyltDQijim+fl7Viq9DN23YAMBLp0/PtxBJklQxA6wqiohWspCnE3g6pbR+N5fwEeBSYAmQyEKsy9j2ksN3AN8FVgODwJ3Aj8hOjnfVF8kubSiZCjx25plnMm3atDF8bPaL+aJFizjjjDNobW0d02dtzz33ZFfTnnZaG69+9Vk1+57dYXf1rJHYs8rZs8rZs8pVs2ddXV1Vqmq3q/oPWkVfAqYBSyJisPjahSmlH470gY0w2rprYIC7N28G4IWTJ4/rUXmO6KycPaucPaucPauco61VCQOsMYqIqWTD8d9MNuy+DQggRcRjZJOkfyeldHuFH72WLGDac9j2PYEnR3pDSulp4PXFydZnA4+TnaAuL9vnYeBlETEZmJZSeiIirijbp/TZewJPDPveu0f53j6gr/S8NP9ra2tr1f6RVs3PGsnvfpetX/7yJlpbG2NquFr3rBHZs8rZs8rZs8pVo2cTrOc784PW2cDbgLcC9wEvAC6OiMdTSpeP8JnjfrT1nS0tDE2ezJ5DQ9xz443cU9Nv2z0c0Vk5e1Y5e1Y5e1a5iTzaWjvPAGsMIuI84EKyE8OfAV8gC416gFnAUcCpwPXFOaQ+lFLaqYvTUkr9EXEH2a+mVxe/r6n4/JIdvLcXWF0cEfZG4MoR9tkCbImImcAryCZuBVhBFmKdRjGwiohpZHcj/NbO1D7eDAzAzTdnj089Nd9aJEkapiY/aAFfAb6UUirNbXlvRMwnC6pGCrDG/WjrW1auhMce48y99uKscX7Ad0Rn5exZ5exZ5exZ5RxtrUoYYI3NicBLyX7hfA1wc0ppU9nrtwHfLU62+m6yMKuS2ZUuAi6PiD8UP+tcYDLZr6hExPeB1SmlC4rPTwb2JQue9gU+TXanyS+XPjAiXkE2QuxB4BCyE9glpc9MKaWIuBj4eEQ8RBZofY7s5PfqCmofN+65BzZvhunT4aij8q5GkjSeVXs6gRr+oNUJDA17yyCj3KG6IUZbb8pO0V42c2bD/MPSEZ2Vs2eVs2eVs2eVc7S1doYB1hiklN5SehwRPwKeB2waYb8+sltYV/r5V0TEHsBnyU6G7wZemVIqzYMxj21PPjvIbp19ELAZWAi8I6W0oWyf6WS/ou4HrAeuIpvzovyC4S+TBWXfAWYAvy1+7/AJYRvCrbdm6xe+EJqbt7+vJEnD1XA6gZKq/6BFNnL8wohYRXYJ4bFkE7h/dxdrrGt9Q0PcVvxl/lQncJckaVwywKqe24ED2XZ4/pillC5hlF9YU0ovH/b8N8CRO/i8KxnhksJh+yTgk8Wl4d12W7Y+qaFvHC5JqoVaTidQUqMftD5ENsL6m8DcYs3/UvyOhrN482b6U2JWSwuHTpqUdzmSJGkXGGBVz9eBL0TEX6SUHs27GO2824u/hxtgSZJ2Qa2nEwBq8oPWJrKRXOdWWst4dHvx8sGTpk179vJHSZI0vhhgVc8VxfV9EfHfwK+Bu4B7U0r9uVWl7erqggceyB6feGK+tUiSxp9aTyeg6ihdPnji1Kk5VyJJknaVAVb1HAgcQ3Yb6mPI7uJzADAQEQ+mlJ6fX2kazR13QEowfz7sOfz+TpIkVaYm0wlo7G4rjcAywJIkadwywKqSlNJKYCXw36VtxUldXwAYXtUp57+SJFWR0wnUoa6BAZZ0dwNw4rRpOVcjSZJ2lQFWDRXnl/if4qI6ZIAlSaoipxOoQ3ds2kQC5rW3s2dbW97lSJKkXWSApQnNCdwlSVXkdAJ16LayCdwlSdL4ZYClCWvdOni0eIHHC16QaymSpAbgdAL16XYncJckqSEYYGnCuueebH3QQeCPspKkWnA6gfzdvXkzAMdNmZJzJZIkaSya8i5AykspwDrmmHzrkCRJtbF5YICHe3sBOMYAS5Kkcc0AazeIiKGI+GVEHJ93LdrKAEuSpMZ275YtAOzd1sYeTuAuSdK4ZoC1e7wHuAn4Rt6FaCsDLEmSGts9xcsHnz95cs6VSJKksTLA2g1SSt9LKX06pfTCvGtRplCA++/PHhtgSZJqzdHY+VhcHIHl5YOSJI1/TuJeBRExh2yU1SnAXsXNTwI3A99LKT2dV20a2ZIl0N8PU6fCAQfkXY0kaQJ4D3AA2Whsf9DaTUojsAywJEka/xyBNUYRcSKwFPgwsJHsUsGbio8/DCyJiBPyq1AjKV0++PznQ0S+tUiSGp+jsXe/oZSeHYHlJYSSJI1/jsAau68D/wGck1JK5S9ERADfLu5zSg61aRSLF2drLx+UJKkxPdLby+bBQdoiOLyzM+9yJEnSGBlgjd0xwLuGh1cAKaUUEV8F7tr9ZWl7SvNfHXVUvnVIkhqD0wnUn8XFywefN3kyrU1edCBJ0njn0XzsngRO2s7rJwFP7aZatJMeeCBbH3FEvnVIksY/pxOoT/d3dwNZgCVJksY/R2CN3T8C3yneVehGtoZVewKnAe8F/ndOtWkEvb2wYkX22ABLklQFTidQhx4oBlhHePmgJEkNwQBrjFJK34iItcDfAn8NNBdfGgTuILu88Mq86tNzLV0KKcHMmTB3bt7VSJIagNMJ1KElBliSJDUUA6wqSCldAVwREa3AnOLmtSmlQo5laRSlywcXLPAOhJKkqihNJ7BklNedTmA3Syk9G2AtMMCSJKkhGGBVUTGweiLvOrR9S4r/vPDyQUlSlTidQJ1Z3dfH5sFBWiI4ZNKkvMuRJElVYIA1BhExL6W0qoL9900pra5lTdoxJ3CXJFWT0wnUn9L8Vwd3dHgHQkmSGoRH9LG5PSL+pXj3oRFFxPSIeG9E/BF4426sTaMojcBasCDfOiRJjSOldEVK6YVAJ7BvcelMKb3Q8Gr3e3b+K+9AKElSw3AE1tgcCVwILIqIXrJfWR8HeoGZxdefB9wJnJ9SWphXocoMDsKDD2aPHYElSao2pxOoD85/JUlS43EE1hiklNYB5wN7A38DPEQ2ifuhxV1+CByfUjrF8Ko+rFwJvb3Q3g4HHJB3NZKk8S4i5lW4/761qkVbPeAdCCVJajiOwBq7lcA/A/+SUvpJ3sVo+x56KFsfcgg0N29/X0mSdsLtEXE18K8ppdtH2iEipgNnAx8BvkN23qAacgSWJEmNxwBr7C4mG3318Yj4LnBxSmlFviVpNMuWZetDDsm3DklSw3A6gTqzaWCAJ/r7ATjcAEuSpIbhJYRjlFL6CnAQ8D7ghcDSiPhJRJycb2UaycMPZ+uDD863DklSY0gprUspnUc2ncAHcTqB3C3v7QVgTmsr01v8rVaSpEbhUb0KUkqDwI+AH0XEqcB5wO8i4lbgH4GrU0opzxqVMcCSJNVCSqknIjamlM7Nu5aJ7uGeHgAO7ujIuRJJklRNjsCqspTS/6SU/hw4jOwygu8BS3MtSs8ywJIk1dA1EXFRRLTlXchEtqwUYE2alHMlkiSpmhyBNUYR8Rlg+ijLDKCT7BJD5SwlWL48e2yAJUmqgZcCPwBOj4i3ppT+OHyHiNgbuCSl9MbdXt0E8bABliRJDckAa+w+QTZR6/fIJmjdCHQVl41la+XsiSegpye7++D8+XlXI0lqNCmlWyPiOLK7DN4eERemlC4CiIgmYAHwt8CpOZbZ8AywJElqTAZYY3ca8FHgPcCPgX8c6RdX5a90+eC8edDamm8tkqTGlFLaHBEfBbqBr0TEW8imbDgSaAdWAhfkWGLDe7g4ibtzYEmS1FicA2uMUkq/Sin9GXAM0AfcGhHXRsRpOZemYZz/SpJUSxHxVxGxClgLvAu4DRgAjgX+FZiVUjowpfT/8quysfUPDbGqFGA5AkuSpIZigFUlKaUHU0rvBw4Afg/8MCLuioi3RURzvtUJDLAkSTX3BeAastFWU1NKp6SUTiEbqf1XwEUR0ZlngY1uZW8vQ0BnUxN7tTmXviRJjcQAq8pSSk+nlD5NNs/Ff5LNg7E816IEwLJl2doAS5JUI78GPl38USuVNqaUvgqcBJwALI6Ik3Oqr+GV5r86aNIkIiLnaiRJUjU5B9YYRcRVjHwHwlagdOY0I5fitI3SCKxDDsm3DklSY0opnb2d1+6NiBOB/wvcRDYflqrM+a8kSWpcBlhj1w08DmzYwaKcrViRrQ86KN86JEkTU0qpDzg3Iq7Ju5ZG5R0IJUlqXAZYY5RSekfeNWjHurth7drs8fz5+dYiSZrYUkqL8q6hUS0vu4RQkiQ1FufA0oSwalW2njoVpk/PtxZJklQbq/r6ADjASwglSWo4BliaEEoB1vz54JyukiQ1plXFObDmtTvFmCRJjcYASxPCypXZet68fOuQJEm1sWVwkHUDAwDMcwSWJEkNxwBLE0L5CCxJktR4SqOvpjc3M73FaV4lSWo0BliaEByBJUlSYyvNf+XoK0mSGpMBVp2LiA9GxCMR0RsRt0bESdvZtzUiPhkRDxf3vyciXjlsn+aI+FxErIiInuK+n4jYOjNURHwvItKw5dpa/p21VhqBZYAlSVJjcv4rSZIam+Or61hEvAm4CDgHuBU4F7guIg5PKa0Z4S2fB94OvBdYArwC+K+IeFFK6a7iPn8PfAB4J3AfcAJwGbAR+Oeyz7oWeHfZ874q/Vm58BJCSZIa28pigDXfEViSJDUkR2DVt/OAS1NKl6WU7icLsrqB94yy/zuAL6SUFqaUlqeUvgUsBD5ats+LgJ+mlK5JKT2SUvoJcD0wfGRXX0rpybLlmar+ZbvR4CA8+mj22BFYkqTxqAYjsh8ZYbR1iohv1P6vqQ0vIZQkqbE5AqtORUQbcDzwxdK2lNJQRNwAnDLK29qB3mHbeoCXlD2/GXhfRByWUloaEccUXz9v2PteHhFrgGeAXwIfTymtG6XW9uJ3l0wFKBQKFAqF7fyVO1Z6/1g+Z/VqGBhopbk5scceA4yxpLpXjZ5NNPascvascvasctXs2Xjue41GZJ8INJe95yhgEfAfNfkjdgMvIZQkqbEZYNWvOWQnlk8N2/4UsGCU91wHnBcRNwEPA6cBb2DbE9QvAdOAJRExWHztwpTSD8v2uRb4T2AFcDDwBeAXEXFKSmlwhO+9APjU8I3XX389nZ2d2/0jd9aiRYt2+b1LlswEXsqsWT1cd92uf854M5aeTVT2rHL2rHL2rHLV6Fl3d3cVKsnNsyOyASLiHODVZCOyvzTC/u8A/iGltLD4/FsRcTrZiOy3A6SUni5/Q0R8jOzc4Tc1+Qt2g5XFEVheQihJUmMywGosHwEuJfu1NZGdiF7Gtpccng28DXgr2RxYLwAujojHU0qXA6SUfly2/70Rsbj4WS8Hbhzhe79I9stwyVTgsTPPPJNp06aN6Q8qFAosWrSIM844g9bW1l36jE2bsvnpDz+8g7POOmtM9YwH1ejZRGPPKmfPKmfPKlfNnnV1dVWpqt2rhiOyh3/H24GLUkpplH3qerT1YEo8Vgyw9m5uHtcj7naGIzorZ88qZ88qZ88q52hrVcIAq36tBQaBPYdt3xN4cqQ3FH9NfX1EdACzgcfJfpldXrbbV4AvlYVU90bEfLJRVJeP8rnLI2ItcAgjBFgppT7KJnkv3dCwtbW1av9IG8tnrV6drQ84oInW1okz7Vs1+z9R2LPK2bPK2bPKVaNn47jntRqRXe71wAzge9upo65HW6+LYGDaNJpT4q4bb2RxVSqqf47orJw9q5w9q5w9q5yjrbUzDLDqVEqpPyLuIDvpvBogIpqKzy/ZwXt7gdUR0Qq8Ebiy7OVOYGjYWwbZzoT+EbEfWSD2RGV/RX0o3YHQCdwlSRPEzozILvf/Ab9IKT2+nc+s69HWv+/qgnvvZf+ODl7zkhEHmjUUR3RWzp5Vzp5Vzp5VztHWqoQBVn27CLg8Iv4A3EY2aetkspNQIuL7wOqU0gXF5ycD+wJ3F9efJgumvlz2mT8DLoyIVWSXEB5LNrfGd4ufMYXsF9aryEZ6HVx8/zKyX3THnVKANX9+vnVIkrQLajUiG4DiKOzTyUZojaruR1sPDADZHQgn0j8aHdFZOXtWOXtWOXtWuQk+2lo7yQCrjqWUroiIPYDPAnuRBVOvTCmVLiOYx7ajqTrI7jx0ELAZWAi8I6W0oWyfDwGfA74JzCU7qf2X4ndAdpL8fOCdZJcTPA5cD3yiePI67qxcma0dgSVJGm9qOCK75N3AGuCaKpa9260qzn81zwncJUlqWAZYdS6ldAmjnKCmlF4+7PlvgCN38HmbyEZynTvK6z1kt9tuGI7AkiSNc7UYkV0Kwt4NXJ5SGtgNf0fNrOrN5qyf396+gz0lSdJ4ZYClhtbVBRs3Zo/33z/fWiRJ2hU1GpEN2aWD8yhOIzCerXQEliRJDc8ASw3t0Uez9YwZMGVKrqVIkrTLqj0iu7jf9UBUo768lUZgzXMEliRJDWvUO89JjaAUYDn/lSRJjcs5sCRJanwGWGpopQDLywclSWpMmwYG2FC8C+H+jsCSJKlhGWCpoZUmcDfAkiSpMT1aHH01o6WFqS3OjiFJUqMywFJDcwSWJEmNrRRgOfpKkqTGZoClhmaAJUlSY3u0OIG7AZYkSY3NAEsNzQBLkqTG9qgTuEuSNCEYYKlhpWSAJUlSo1vlJYSSJE0IBlhqWOvWQfGqAvbbL99aJElSbXgJoSRJE4MBlhpWafTV3LngOa0kSY3JSdwlSZoYDLDUsLx8UJKkxpZS2hpgOQeWJEkNzQBLDcsAS5KkxrauUKBnaAiA/RyBJUlSQzPAUsMywJIkqbGVRl/t2dpKe5OntZIkNTKP9GpYBliSJDU2Lx+UJGniMMBSwzLAkiSpsTmBuyRJE4cBlhqWAZYkSY1tVW8vAPMMsCRJangGWGpIQ0OwenX22ABLkqTG5CWEkiRNHAZYakhPPQWFAjQ1wT775F2NJEmqBS8hlCRp4jDAUkNauTJb77MPtLTkW4skSaqNR4uXEBpgSZLU+Ayw1JBWrMjWBx6Ybx2SJKk2BoaGeKw4AmuelxBKktTwDLDUkB55JFsfcECeVUiSpFp5tK+PQaA9gr3b2vIuR5Ik1ZgBlhqSI7AkSWpsK4qXDx7Q0UFTRM7VSJKkWjPAUkMywJIkqbGVAqwDJ03KuRJJkrQ7GGCpIRlgSZLU2Fb09ABwoPNfSZI0IRhgqeEMDsKqVdljAyxJkhrT8tIILAMsSZImBAMsNZzHH4dCAVpbYd99865GkiTVwgoDLEmSJhQDLDWc0uWD8+ZBc3O+tUiSpNp49hJC58CSJGlCMMBSw3n44Wzt5YOSJDWm7sFBnioUADjIEViSJE0IBlhqOA8+mK0PPzzfOiRJUm08XBx9NaOlhZmtrTlXI0mSdgcDLDUcAyxJkhrbku5uABZ0duZciSRJ2l0MsNRwlizJ1gsW5FuHJEmqDQMsSZImHgMsNZRCYescWI7AkiSpMRlgSZI08RhgqaGsWJGFWJ2dsN9+eVcjSZJqoRRgHWGAJUnShGGApYZSPv9Vk/91S5LUcIZScgSWJEkTkP/EV0NZvDhbH3FEvnVIkqTaeLSvj+6hIVoiOLCjI+9yJEnSbmKApYZy113Z+thj861DkiTVxt2bNwPwvM5OWh1uLUnShOFRXw2lFGAdd1y+dUiSpNq4c9MmAI6bOjXnSiRJ0u5kgKWGsXEjLF+ePX7BC3ItRZIk1cidxRFYx02ZknMlkiRpdzLAUsO4++5sPX8+zJqVaymSJKlGHIElSdLEZIClhnHnndna+a8kSWpMT/b18Xh/PwEc4wgsSZImFAMsNQwncJckqbHdVbx8cEFnJ5Obm3OuRpIk7U4GWGoYjsCSJKmxOf+VJEkTlwGWGsKWLfDAA9njE07ItxZJklQbdxTnvzrW+a8kSZpwDLDqXER8MCIeiYjeiLg1Ik7azr6tEfHJiHi4uP89EfHKYfs0R8TnImJFRPQU9/1ERETZPhERn42IJ4r73BARh9by7xyru++GoSHYe+9skSRJjacUYB3vCCxJkiYcA6w6FhFvAi4CPgMcB9wDXBcRc0d5y+eB9wMfAo4Evg38V0SUX1T398AHgL8Bjig+P7/4npLzgQ8D5wAnA1uK39tRnb+s+u64I1sff3y+dUiSpNpY29/Pqr4+wDsQSpI0ERlg1bfzgEtTSpellO4nC5S6gfeMsv87gC+klBamlJanlL4FLAQ+WrbPi4CfppSuSSk9klL6CXA9cBJko6+Ac4HPp5R+mlJaDPwvYB/g9VX/C6ukFGB5+aAkqRFVe0R2cb99I+LfImJdccT1vRFRt0fSO4rzXx02aRLTWlpyrkaSJO1uHv3rVES0AccDXyxtSykNRcQNwCmjvK0d6B22rQd4Sdnzm4H3RcRhKaWlEXFM8fXziq8fCOwF3FD2vRsj4tbi9/54hFrbi99dMhWgUChQKBR29KduV+n9O/qc229vAYJjjhmgUEhj+s7xbmd7pq3sWeXsWeXsWeWq2bPx3PeyEdnnALeS/dB0XUQcnlJaM8JbPg+8HXgvsAR4BdmI7BellO4qfuZM4HfAr4BXAU8DhwLP1Pav2XXPXj7o6CtJkiYkA6z6NQdoBp4atv0pYMEo77kOOC8ibgIeBk4D3lD8nJIvAdOAJRExWHztwpTSD4uv71X2PcO/dy9GdgHwqeEbr7/+ejo7O0d5S2UWLVo06mu9vc08+OCrAXjmmRtYuLCvKt853m2vZxqZPaucPaucPatcNXrW3d1dhUpy8+yIbICIOAd4NdmI7C+NsP87gH9IKS0sPv9WRJxONiL77cVtfw88mlJ6d9n7VtSi+GoxwJIkaWIzwGosHwEuJfu1NZGFWJex7SWHZwNvA94K3Ae8ALg4Ih5PKV2+i9/7RbJfhkumAo+deeaZTJs2bRc/MlMoFFi0aBFnnHEGra2tI+6TTeAezJ6dePvbTxvT9zWCnemZtmXPKmfPKmfPKlfNnnV1dVWpqt2rhiOyX0s2ius/gJcBq4FvppQurVbt1XZn8RJCJ3CXJGliMsCqX2uBQWDPYdv3BJ4c6Q0ppaeB1xcnW58NPE72y+zyst2+AnwppVS6FPDeiJhPNorq8rLP3hN4Ytj33j3K9/YBzw57Kt3QsLW1tWr/SNveZ61cma0PPTT8R2GZavZ/orBnlbNnlbNnlatGz8Zxz2s1Ivsgspu6XAR8ATgR+OeI6B/pB628pwvoHRpiZW+WyR3W3j6uLwmtBi9Jrpw9q5w9q5w9q5zTBagSBlh1KqXUHxF3kJ10Xg0QEU3F55fs4L29wOqIaAXeCFxZ9nInMDTsLYNsndB/BVmIdRrFwCoippHdjfBbu/wH1dBDD2XrQw/Ntw5JkurEzozIbgL+kFL6P8Xnd0XEUWTzbI00IjvX6QIebWoiTZ1KZ0rctmgRUZVvHP+8JLly9qxy9qxy9qxyThegnWGAVd8uAi6PiD8At5FN2jqZ7CSUiPg+sDqldEHx+cnAvmTB077Ap8lOUL9c9pk/Ay6MiFVklxAeSza3xncBUkopIi4GPh4RD5EFWp8jG811da3+0LFYtixbH3JIvnVIklQDtRqR/QRw/7C3PkD2w9dIcp0u4Gfr1sGSJSyYMoVXv+QlI+4zkXhJcuXsWeXsWeXsWeWcLkCVMMCqYymlKyJiD+CzZBOo3w28MqVUuoxgHtuOpuogu/PQQcBmYCHwjpTShrJ9PkQWSH0TmEt2Uvsvxe8o+TJZUPYdYAbw2+L3Dp9Poy4YYEmSGlUNR2T/Djh82FsOA1aO8lm5ThewonhZyKGdnf6jsIyXJFfOnlXOnlXOnlVugk8XoJ1kgFXnUkqXMMoJakrp5cOe/wY4cgeft4lsJNe529knAZ8sLnXPSwglSQ2uFiOyvwrcHBH/hyzYOgl4X3GpOw8VLws5tEqXK0qSpPHHAEvj2pYt8ERxqnlHYEmSGlEtRmSnlG6PiD8nuzTwk2RTBpybUvphbf+aXbOspweAQyZNyrkSSZKUFwMsjWulywdnzYKZM/OtRZKkWqn2iOzifj8Hfl6N+mrtoWKAdagBliRJE1bTjneR6lcpwPLyQUmSGlPv4CCP9mXTbzkCS5KkicsAS+OaE7hLktTYlvf2koBpzc3s4QS9kiRNWAZYGtdKE7gbYEmS1JjK578q3f1QkiRNPAZYGte8hFCSpMbm/FeSJAkMsDTOeQmhJEmN7aHubsD5ryRJmugMsDRudXfD6tXZYwMsSZIaU+kSwkM7O3OuRJIk5ckAS+PWww9n65kzYfbsfGuRJEm18VDZHFiSJGniMsDSuOXlg5IkNbbewUEe7esDnANLkqSJzgBL41bpDoRO4C5JUmNa3ttLAqY1N7NHa2ve5UiSpBwZYGnccgSWJEmNbVnZ5YMRkXM1kiQpTwZYGrdKI7AMsCRJakyl+a+8fFCSJBlgadwqjcDyEkJJkhrTQ93dgBO4S5IkAyyNUz098Nhj2WNHYEmS1JhKlxAe2tmZcyWSJClvBlgalx5+OFtPnw6zZ+dbiyRJqo2HyubAkiRJE5sBlsal8gncndNVkqTG0zs4yKN9fYABliRJMsDSOPXgg9n68MPzrUOSJNXGsp4eEjC9uZm5ra15lyNJknJmgKVxaenSbH3YYfnWIUmSauPB4uWDh3d2Eg63liRpwjPA0rjkCCxJkhrbg8U7EB7mBO6SJAkDLI1TpRFYBliSJDWmUoB1uPNfSZIkDLA0Dj3zDDz9dPb40EPzrUWSJNXG0rJLCCVJkgywNO788Y/Zev/9YcqUfGuRJEnVN5QSf9yyBYAjDLAkSRIGWBqH7rwzWx93XL51SJKk2ljW08PmwUE6mppYYIAlSZIwwNI4dMcd2doAS5KkxnTnpk0AHDN5Mi1Nnq5KkiQDLI1Dt92WrQ2wJElqTLcVA6zjpk7NuRJJklQvDLA0rjz2GDz4IDQ1wUteknc1kiSpFhY98wwAL58xI99CJElS3TDA0rjy859n6xNPBM9pJUlqPCt6evjjli0E8Kce7CVJUpEBlsaNwUH45jezx2efnW8tkiSpNi5ZvRqAM2bOZE5bW87VSJKketGSdwHS9lxySRM33XQk117bxH33wb33wvTp8M535l2ZJEmqhl+sX88P2ttZtHw5Tw8M8JOnnwbgI/vtl3NlkiSpnhhgqa5demkTDzxw6LPPm5vh0kth9uwci5IkSVVz44YNXNXRAU888ey2c/bZh7M82EuSpDIGWKprb3vbELfdtoKjjz6QmTObefWrYcGCvKuSJEnVcur06Sx/5BGOOvhgpre2csLUqZw2c2beZUmSpDpjgKW6dv75QyxceB9nnTWf1tbmvMuRJElV9rrZs2nt7eWs+fNpbW3NuxxJklSnnMRdkiRJkiRJdc0AS5IkSZIkSXXNAEuSJEmSJEl1zQBLkiRJkiRJdc0AS5IkSZIkSXXNAEuSJEmSJEl1zQBLkiRJkiRJdc0AS5IkSZIkSXXNAEuSJEmSJEl1rSXvAtS4urq6xvwZhUKB7u5uurq6aG1trUJVjc+eVc6eVc6eVc6eVa6aPavGMUnP5bE+H/ascvascvascvasch7rVYlIKeVdgxpMROwLPJZ3HZIkjWC/lNLqvIsY7zzWS5LqmMf6BmWApaqLiAD2ATZV4eOmkp0g71elz5sI7Fnl7Fnl7Fnl7Fnlqt2zqcDjyZOfMfNYnzt7Vjl7Vjl7Vjl7VjmP9dppXkKoqiv+j0VVEu/s/BiATSklx4TuBHtWOXtWOXtWOXtWuRr0zL5Xicf6fNmzytmzytmzytmzynmsVyWcxF2SJEmSJEl1zQBLkiRJkiRJdc0AS/WuD/hMca2dY88qZ88qZ88qZ88qZ88mBv/vXDl7Vjl7Vjl7Vjl7Vjl7pp3mJO6SJEmSJEmqa47AkiRJkiRJUl0zwJIkSZIkSVJdM8CSJEmSJElSXTPAkiRJkiRJUl0zwFJdi4gPRsQjEdEbEbdGxEl511QvIuKlEfGziHg8IlJEvH7Y6xERn42IJyKiJyJuiIhDcyo3dxFxQUTcHhGbImJNRFwdEYcP26cjIr4REesiYnNEXBURe+ZVc94i4gMRsTgiuorLLRHxqrLX7dcORMTHiv//eXHZNvtWJiI+XexR+bKk7HX71eA81o/OY31lPNZXzmP92Hms3zGP9aoWAyzVrYh4E3AR2W1VjwPuAa6LiLm5FlY/JpP15IOjvH4+8GHgHOBkYAtZ/zp2T3l152XAN4AXAmcArcD1ETG5bJ+vAq8B/rK4/z7Af+7mOuvJY8DHgOOBE4BfAj+NiOcVX7df2xERJwLvBxYPe8m+Pdd9wN5ly0vKXrNfDcxj/Q55rK+Mx/rKeawfA4/1FfFYr7FLKbm41OUC3ApcUva8CVgNfCzv2uptARLw+rLnATwB/O+ybdOBXuDNeddbDwuwR7FvLy3rTz/wF2X7LCju88K8662XBVgP/H/2a4d9mgIsBU4Hfg1cXNxu357bq08Dd4/ymv1q8MVjfUW98lhfec881u9a3zzW71yfPNbvfK881rtUZXEElupSRLSR/RJ0Q2lbSmmo+PyUvOoaRw4E9mLb/m0k+4eC/ctML67XF9fHk/1SW96zJcAq7BkR0RwRbyYbDXAL9mtHvgFck1K6Ydh2+zayQ4uXSC2PiB9GxLzidvvVwDzWj5nH+h3zWF8Bj/UV81hfGY/1GrOWvAuQRjEHaAaeGrb9KbJEXtu3V3E9Uv/2YoKLiCbgYuB3KaU/FjfvBfSnlDYM231C9ywijiY7ie0ANgN/nlK6PyJegP0aUfHk/zjgxBFe9r+z57oVeBfwINklBZ8C/icijsJ+NTqP9WPjsX47PNbvPI/1lfNYXzGP9aoKAyxJE9E3gKPY9tp7jexB4AVkv2L/BXB5RLws14rqWETsD3wNOCOl1Jt3PeNBSukXZU8XR8StwErgbKAnn6okNQCP9TvPY30FPNZXzmO9qsVLCFWv1gKDwPC7T+wJPLn7yxl3Sj2yf8NExCXAnwF/klJ6rOylJ4G2iJgx7C0Tumcppf6U0rKU0h0ppQvIJhP+CPZrNMcDc4E7I2IgIgbIJiP9cPHxU9i37Sr+ArsUOAT/O2t0HuvHxmP9KDzWV8ZjfcU81o+Rx3rtKgMs1aWUUj9wB3BaaVtxKPhpZEOctX0ryP4Hv7x/08juUDQh+1e81fglwJ8Df5pSWjFslzuAAtv27HBgHhO0Z6NoAtqxX6O5ETia7Jfs0vIH4Idlj+3bdkTEFOBgssmp/e+sgXmsHzOP9cN4rK8aj/Xb57F+jDzWa1d5CaHq2UVkQ5j/ANwGnEs2qeRleRZVL4r/w39I2aYDi3MVrE8prYqIi4GPR8RDZCe5nwMeB67ezaXWi28AbwVeB2yKiNI19RtTSj0ppY0R8f+AiyJiPdAFfB24JaX0+3xKzldEfBH4BdkkmlPJ+vdy4BX2a2QppU3AH8u3RcQWYF1pDhb7tq2I+EfgZ2SXEuwDfIZsVM6P/O9sQvBYvx0e6yvmsb5CHusr57G+ch7rVS0GWKpbKaUrImIP4LNkE/jdDbwypTR8stKJ6gTgV2XPLyquLyebJPHLZP8I+A4wA/gtWf8m6rX6Hyiufz1s+7uB7xUf/y0wBFxF9svjdcBf74ba6tVc4Ptkk21uBBaTndAuKr5uv3aNfdvWfsCPgNnA02T/W/XClNLTxdftVwPzWL9DHusr47G+ch7ra8O+bctjvaoiUkp51yBJkiRJkiSNyjmwJEmSJEmSVNcMsCRJkiRJklTXDLAkSZIkSZJU1wywJEmSJEmSVNcMsCRJkiRJklTXDLAkSZIkSZJU1wywJEmSJEmSVNcMsCRJkiRJklTXDLAkSZIkSZJU1wywJGmciIiWvGuQJEm147FekkZngCVJdSgiDoiIFBFnR8T/REQf8Nq865IkSdXhsV6SKmPCL0n16Zji+u+A/wOsAJ7OrxxJklRlHuslqQIGWJJUn14AbAH+MqX0SL6lSJKkGngBHuslaad5CaEk1adjgP/2hFaSpIblsV6SKmCAJUn16QXAr3OuQZIk1c4L8FgvSTvNAEuS6kxETAMOAO7KuRRJklQDHuslqXIGWJJUf44BBoF78y5EkiTVhMd6SaqQAZYk1Z9jgAdTSr15FyJJkmrCY70kVShSSnnXIEmSJEmSJI3KEViSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmsGWJIkSZIkSaprBliSJEmSJEmqawZYkiRJkiRJqmv/P8D1jDdZ+VD6AAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"# Compute phi\n",
"phi = np.log(psi)\n",
"\n",
"# Compute W\n",
"W = psi**(-2)\n",
"\n",
"# Compute chi\n",
"chi = psi**(-4)\n",
"\n",
"f = plt.figure(figsize=(12,8),dpi=100)\n",
"\n",
"ax = f.add_subplot(221)\n",
"ax.set_title(r\"Conformal factor $\\psi(0,r)$\")\n",
"ax.set_ylabel(r\"$\\psi(0,r)$\")\n",
"ax.plot(r,psi,'k-')\n",
"ax.grid()\n",
"\n",
"ax2 = f.add_subplot(222)\n",
"ax2.set_title(r\"Evolved conformal factor $\\phi(0,r)$\")\n",
"ax2.set_ylabel(r\"$\\phi(0,r)$\")\n",
"ax2.plot(r,phi,'r-')\n",
"ax2.grid()\n",
"\n",
"ax3 = f.add_subplot(223)\n",
"ax3.set_title(r\"Evolved conformal factor $W(0,r)$\")\n",
"ax3.set_xlabel(r\"$r$\")\n",
"ax3.set_ylabel(r\"$W(0,r)$\")\n",
"ax3.plot(r,W,'b-')\n",
"ax3.grid()\n",
"\n",
"ax4 = f.add_subplot(224)\n",
"ax4.set_title(r\"Evolved conformal factor $\\chi(0,r)$\")\n",
"ax4.set_xlabel(r\"$r$\")\n",
"ax4.set_ylabel(r\"$\\chi(0,r)$\")\n",
"ax4.plot(r,chi,'c-')\n",
"ax4.grid()\n",
"\n",
"outfile = os.path.join(Ccodesdir,\"cfs_scalarfield_id.png\")\n",
"plt.savefig(outfile)\n",
"plt.close(f)\n",
"\n",
"# Display the figure\n",
"from IPython.display import Image\n",
"Image(outfile)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 1.d The lapse function: $\\tilde\\alpha$ \\[Back to [top](#toc)\\]\n",
"$$\\label{id_lapse_function}$$\n",
"\n",
"There are two common initial conditions for $\\tilde\\alpha$. The first one is eq. (44) of [A&C](https://arxiv.org/pdf/1508.01614.pdf), namely setting the lapse to unity\n",
"\n",
"$$\n",
"\\tilde\\alpha = 1\\ .\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:10.020202Z",
"iopub.status.busy": "2021-06-15T10:05:10.019698Z",
"iopub.status.idle": "2021-06-15T10:05:10.021704Z",
"shell.execute_reply": "2021-06-15T10:05:10.021195Z"
}
},
"outputs": [],
"source": [
"# Set the unity lapse initial condition\n",
"alpha_unity = np.ones(NR)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The second one is discussed in the last paragraph of section II.B in [Baumgarte](https://arxiv.org/pdf/1807.10342.pdf), which is to set the \"pre-collapsed lapse\"\n",
"\n",
"$$\n",
"\\tilde\\alpha = \\psi^{-2}\\ .\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:10.024115Z",
"iopub.status.busy": "2021-06-15T10:05:10.023725Z",
"iopub.status.idle": "2021-06-15T10:05:10.026733Z",
"shell.execute_reply": "2021-06-15T10:05:10.026335Z"
}
},
"outputs": [],
"source": [
"# Set the \"pre-collapsed lapse\" initial condition\n",
"alpha_precollapsed = psi**(-2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 1.e Outputting the initial data to file \\[Back to [top](#toc)\\]\n",
"$$\\label{id_output}$$"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:10.031335Z",
"iopub.status.busy": "2021-06-15T10:05:10.030923Z",
"iopub.status.idle": "2021-06-15T10:05:10.245788Z",
"shell.execute_reply": "2021-06-15T10:05:10.245277Z"
}
},
"outputs": [],
"source": [
"# Check to see which version of Python is being used\n",
"# For a machine running the final release of Python 3.7.1,\n",
"# sys.version_info should return the tuple [3,7,1,'final',0]\n",
"if sys.version_info[0] == 3:\n",
" np.savetxt(os.path.join(Ccodesdir,\"outputSFID_unity_lapse.txt\"), list(zip( r, ID_sf, psi**4, alpha_unity )),\n",
" fmt=\"%.15e\")\n",
" np.savetxt(os.path.join(Ccodesdir,\"outputSFID_precollapsed_lapse.txt\"), list(zip( r, ID_sf, psi**4, alpha_precollapsed )),\n",
" fmt=\"%.15e\")\n",
"elif sys.version_info[0] == 2:\n",
" np.savetxt(os.path.join(Ccodesdir,\"outputSFID_unity_lapse.txt\"), zip( r, ID_sf, psi**4, alpha_unity ),\n",
" fmt=\"%.15e\")\n",
" np.savetxt(os.path.join(Ccodesdir,\"outputSFID_precollapsed_lapse.txt\"), zip( r, ID_sf, psi**4, alpha_precollapsed ),\n",
" fmt=\"%.15e\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 2: Interpolating the initial data file as needed \\[Back to [top](#toc)\\]\n",
"$$\\label{id_interpolation_files}$$\n",
"\n",
"In order to use the initial data file properly, we must tell the program how to interpolate the values we just computed to the values of $r$ in our numerical grid. We do this by creating two C functions: one that interpolates the ADM quantities, $\\left\\{\\gamma_{ij},K_{ij},\\alpha,\\beta^{i},B^{i}\\right\\}$, and one that interpolates the scalar field quantities, $\\left\\{\\varphi,\\Pi\\right\\}$. The two files written below use the scalarfield_interpolate_1D( ) function, which is defined in the [ScalarField/ScalarField_interp.h](../edit/ScalarField/ScalarField_interp.h) file. This function performs a Lagrange polynomial interpolation between the initial data file and the numerical grid used during the simulation."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:10.249253Z",
"iopub.status.busy": "2021-06-15T10:05:10.248731Z",
"iopub.status.idle": "2021-06-15T10:05:10.250854Z",
"shell.execute_reply": "2021-06-15T10:05:10.250410Z"
}
},
"outputs": [],
"source": [
"def ID_scalarfield_ADM_quantities(Ccodesdir=\".\",new_way=False):\n",
" includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
" desc = \"\"\"(c) 2021 Leo Werneck\n",
"This function takes as input either (x,y,z) or (r,th,ph) and outputs\n",
"all ADM quantities in the Cartesian or Spherical basis, respectively.\n",
"\"\"\"\n",
" c_type = \"void\"\n",
" name = \"ID_scalarfield_ADM_quantities\"\n",
" params = \"\"\"const REAL xyz_or_rthph[3],const ID_inputs other_inputs,\n",
" REAL *restrict gammaDD00,REAL *restrict gammaDD01,REAL *restrict gammaDD02,\n",
" REAL *restrict gammaDD11,REAL *restrict gammaDD12,REAL *restrict gammaDD22,\n",
" REAL *restrict KDD00,REAL *restrict KDD01,REAL *restrict KDD02,\n",
" REAL *restrict KDD11,REAL *restrict KDD12,REAL *restrict KDD22,\n",
" REAL *restrict alpha,\n",
" REAL *restrict betaU0,REAL *restrict betaU1,REAL *restrict betaU2,\n",
" REAL *restrict BU0,REAL *restrict BU1,REAL *restrict BU2\"\"\"\n",
" body = \"\"\"\n",
" const REAL r = xyz_or_rthph[0];\n",
" const REAL th = xyz_or_rthph[1];\n",
" const REAL ph = xyz_or_rthph[2];\n",
"\n",
" REAL sf_star,psi4_star,alpha_star;\n",
"\n",
" scalarfield_interpolate_1D(r,\n",
" other_inputs.interp_stencil_size,\n",
" other_inputs.numlines_in_file,\n",
" other_inputs.r_arr,\n",
" other_inputs.sf_arr,\n",
" other_inputs.psi4_arr,\n",
" other_inputs.alpha_arr,\n",
" &sf_star,&psi4_star,&alpha_star);\n",
"\n",
" // Update alpha\n",
" *alpha = alpha_star;\n",
" // gamma_{rr} = psi^4\n",
" *gammaDD00 = psi4_star;\n",
" // gamma_{thth} = psi^4 r^2\n",
" *gammaDD11 = psi4_star*r*r;\n",
" // gamma_{phph} = psi^4 r^2 sin^2(th)\n",
" *gammaDD22 = psi4_star*r*r*sin(th)*sin(th);\n",
"\n",
" // All other quantities ARE ZERO:\n",
" *gammaDD01 = 0.0; *gammaDD02 = 0.0;\n",
" /**/ *gammaDD12 = 0.0;\n",
"\n",
" *KDD00 = 0.0; *KDD01 = 0.0; *KDD02 = 0.0;\n",
" /**/ *KDD11 = 0.0; *KDD12 = 0.0;\n",
" /**/ *KDD22 = 0.0;\n",
"\n",
" *betaU0 = 0.0; *betaU1 = 0.0; *betaU2 = 0.0;\n",
"\n",
" *BU0 = 0.0; *BU1 = 0.0; *BU2 = 0.0;\n",
"\"\"\"\n",
" if new_way == True:\n",
" outC.add_to_Cfunction_dict(includes=includes,desc=desc,c_type=c_type,name=name,\n",
" params=params,body=body,enableCparameters=False)\n",
" else:\n",
" outfile = os.path.join(Ccodesdir,\"ID_scalarfield_ADM_quantities-validation.h\")\n",
" outC.outCfunction(outfile=outfile,\n",
" includes=None,desc=desc,c_type=c_type,name=name,\n",
" params=params,body=body,enableCparameters=False)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def ID_scalarfield_spherical(Ccodesdir=\".\",new_way=False):\n",
" includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
" desc = \"\"\"(c) 2021 Leo Werneck\n",
"This function takes as input either (x,y,z) or (r,th,ph) and outputs all\n",
"scalar field quantities in the Cartesian or Spherical basis, respectively.\n",
"\"\"\"\n",
" c_type = \"void\"\n",
" name = \"ID_scalarfield_spherical\"\n",
" params = \"const REAL xyz_or_rthph[3],const ID_inputs other_inputs,REAL *restrict sf,REAL *restrict sfM\"\n",
" body = \"\"\"\n",
" const REAL r = xyz_or_rthph[0];\n",
" const REAL th = xyz_or_rthph[1];\n",
" const REAL ph = xyz_or_rthph[2];\n",
"\n",
" REAL sf_star,psi4_star,alpha_star;\n",
"\n",
" scalarfield_interpolate_1D(r,\n",
" other_inputs.interp_stencil_size,\n",
" other_inputs.numlines_in_file,\n",
" other_inputs.r_arr,\n",
" other_inputs.sf_arr,\n",
" other_inputs.psi4_arr,\n",
" other_inputs.alpha_arr,\n",
" &sf_star,&psi4_star,&alpha_star);\n",
"\n",
" // Update varphi\n",
" *sf = sf_star;\n",
" // Update Pi\n",
" *sfM = 0;\n",
"\"\"\"\n",
" if new_way == True:\n",
" outC.add_to_Cfunction_dict(includes=includes,desc=desc,c_type=c_type,name=name,\n",
" params=params,body=body,enableCparameters=False)\n",
" else:\n",
" outfile = os.path.join(Ccodesdir,\"ID_scalarfield_spherical-validation.h\")\n",
" outC.outCfunction(outfile=outfile,\n",
" includes=None,desc=desc,c_type=c_type,name=name,\n",
" params=params,body=body,enableCparameters=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 3: Converting Spherical initial data to Curvilinear initial data \\[Back to [top](#toc)\\]\n",
"$$\\label{id_sph_to_curvilinear}$$\n",
"\n",
"In this tutorial module we have explained how to obtain spherically symmetric, time-symmetric initial data for the collapse of a massless scalar field in Spherical coordinates (see [Step 1](#initial_data)). We have also explained how to interpolate the initial data file to the numerical grid we will use during the simulation (see [Step 2](#id_interpolation_files)).\n",
"\n",
"NRPy+ is capable of generating the BSSN evolution equations in many different Curvilinear coordinates (for example SinhSpherical coordinates, which are of particular interest for this problem). Therefore, it is essential that we convert the Spherical initial data generated here to any Curvilinear system supported by NRPy+.\n",
"\n",
"We start by calling the reference_metric() function within the [reference_metric.py](../edit/reference_metric.py) NRPy+ module. This will set up a variety of useful quantities for us."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then the code below interpolate the values of the Spherical grid $\\left\\{r,\\theta,\\phi\\right\\}$ to the Curvilinear grid $\\left\\{{\\rm xx0,xx1,xx2}\\right\\}$."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:10.499338Z",
"iopub.status.busy": "2021-06-15T10:05:10.498812Z",
"iopub.status.idle": "2021-06-15T10:05:10.506013Z",
"shell.execute_reply": "2021-06-15T10:05:10.506415Z"
}
},
"outputs": [],
"source": [
"def ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2(Ccodesdir=\".\",pointer_to_ID_inputs=False,new_way=False):\n",
"\n",
" rfm.reference_metric()\n",
"\n",
" rthph = outC.outputC(rfm.xxSph[0:3],[\"rthph[0]\", \"rthph[1]\", \"rthph[2]\"],\n",
" \"returnstring\", \"includebraces=False,outCverbose=False,preindent=1\")\n",
"\n",
" includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
" desc = \"\"\"(c) 2021 Leo Werneck\n",
"This function takes as input either (x,y,z) or (r,th,ph) and outputs all\n",
"scalar field quantities in the Cartesian or Spherical basis, respectively.\n",
"\"\"\"\n",
" c_type = \"void\"\n",
" name = \"ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2\"\n",
" params = \"const paramstruct *restrict params,const REAL xx0xx1xx2[3],\\n\"\n",
" if pointer_to_ID_inputs == True:\n",
" params += \"ID_inputs *other_inputs,\\n\"\n",
" else:\n",
" params += \"ID_inputs other_inputs,\\n\"\n",
" params += \"REAL *restrict sf, REAL *restrict sfM\"\n",
" body = \"\"\"\n",
" const REAL xx0 = xx0xx1xx2[0];\n",
" const REAL xx1 = xx0xx1xx2[1];\n",
" const REAL xx2 = xx0xx1xx2[2];\n",
" REAL rthph[3];\n",
"\"\"\"+rthph+\"\"\"\n",
" ID_scalarfield_spherical(rthph,other_inputs,sf,sfM);\n",
"\"\"\"\n",
"\n",
" if new_way == True:\n",
" outC.add_to_Cfunction_dict(includes=includes,desc=desc,c_type=c_type,name=name,\n",
" params=params,body=body)\n",
" else:\n",
" outfile = os.path.join(Ccodesdir,\"ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2-validation.h\")\n",
" outC.outCfunction(outfile=outfile,\n",
" includes=None,desc=desc,c_type=c_type,name=name,\n",
" params=params,body=body)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we create the driver function which puts everything together using OpenMP."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:10.509545Z",
"iopub.status.busy": "2021-06-15T10:05:10.509031Z",
"iopub.status.idle": "2021-06-15T10:05:10.511088Z",
"shell.execute_reply": "2021-06-15T10:05:10.510624Z"
}
},
"outputs": [],
"source": [
"def ID_scalarfield(Ccodesdir=\".\",new_way=False):\n",
" includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
" desc = \"\"\"(c) 2021 Leo Werneck\n",
"This is the scalar field initial data driver functiono.\n",
"\"\"\"\n",
" c_type = \"void\"\n",
" name = \"ID_scalarfield\"\n",
" params = \"\"\"const paramstruct *restrict params,REAL *restrict xx[3],\n",
" ID_inputs other_inputs,REAL *restrict in_gfs\"\"\"\n",
" body = \"\"\"\n",
" const int idx = IDX3S(i0,i1,i2);\n",
" const REAL xx0xx1xx2[3] = {xx0,xx1,xx2};\n",
" ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2(params,xx0xx1xx2,other_inputs,\n",
" &in_gfs[IDX4ptS(SFGF,idx)],\n",
" &in_gfs[IDX4ptS(SFMGF,idx)]);\n",
"\"\"\"\n",
" loopopts = \"AllPoints,Read_xxs\"\n",
" if new_way == True:\n",
" outC.add_to_Cfunction_dict(includes=includes,desc=desc,c_type=c_type,name=name,\n",
" params=params,body=body,loopopts=loopopts)\n",
" else:\n",
" outfile = os.path.join(Ccodesdir,\"ID_scalarfield-validation.h\")\n",
" outC.outCfunction(outfile=outfile,\n",
" includes=None,desc=desc,c_type=c_type,name=name,\n",
" params=params,body=body,loopopts=loopopts)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(Ccodesdir=\".\",pointer_to_ID_inputs=False,new_way=False):\n",
" ID_scalarfield_ADM_quantities(Ccodesdir=Ccodesdir,new_way=new_way)\n",
" ID_scalarfield_spherical(Ccodesdir=Ccodesdir,new_way=new_way)\n",
" ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2(Ccodesdir=Ccodesdir,pointer_to_ID_inputs=pointer_to_ID_inputs,new_way=new_way)\n",
" ID_scalarfield(Ccodesdir=Ccodesdir,new_way=new_way)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Output C function ID_scalarfield_ADM_quantities() to file ScalarFieldID_validation/ID_scalarfield_ADM_quantities-validation.h\n",
"Output C function ID_scalarfield_spherical() to file ScalarFieldID_validation/ID_scalarfield_spherical-validation.h\n",
"Output C function ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2() to file ScalarFieldID_validation/ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2-validation.h\n",
"Output C function ID_scalarfield() to file ScalarFieldID_validation/ID_scalarfield-validation.h\n"
]
}
],
"source": [
"NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(Ccodesdir=Ccodesdir)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 4: Validation of this tutorial against the [ScalarField/ScalarField_InitialData.py](../edit/ScalarField/ScalarField_InitialData.py) module \\[Back to [top](#toc)\\]\n",
"$$\\label{validation}$$\n",
"\n",
"First we load the [ScalarField/ScalarField_InitialData.py](../edit/ScalarField/ScalarField_InitialData.py) module and compute everything by using the scalarfield_initial_data( ) function, which should do exactly the same as we have done in this tutorial."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:10.515012Z",
"iopub.status.busy": "2021-06-15T10:05:10.514505Z",
"iopub.status.idle": "2021-06-15T10:05:10.849742Z",
"shell.execute_reply": "2021-06-15T10:05:10.850162Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Generated the ADM initial data for the gravitational collapse \n",
"of a massless scalar field in Spherical coordinates.\n",
"\n",
"Type of initial condition: Scalar field: \"Gaussian\" Shell\n",
" ADM quantities: Time-symmetric\n",
" Lapse condition: Unity\n",
"Parameters: amplitude = 0.1,\n",
" center = 0,\n",
" width = 1,\n",
" domain size = 50,\n",
" number of points = 30000,\n",
" Initial data file = ScalarFieldID_validation/outputSFID_unity_lapse-validation.txt.\n",
"\n",
"Generated the ADM initial data for the gravitational collapse \n",
"of a massless scalar field in Spherical coordinates.\n",
"\n",
"Type of initial condition: Scalar field: \"Gaussian\" Shell\n",
" ADM quantities: Time-symmetric\n",
" Lapse condition: Pre-collapsed\n",
"Parameters: amplitude = 0.1,\n",
" center = 0,\n",
" width = 1,\n",
" domain size = 50,\n",
" number of points = 30000,\n",
" Initial data file = ScalarFieldID_validation/outputSFID_precollapsed_lapse-validation.txt.\n",
"\n",
"Output C function ID_scalarfield_ADM_quantities() to file ScalarFieldID_validation/ID_scalarfield_ADM_quantities.h\n",
"Output C function ID_scalarfield_spherical() to file ScalarFieldID_validation/ID_scalarfield_spherical.h\n",
"Output C function ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2() to file ScalarFieldID_validation/ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2.h\n",
"Output C function ID_scalarfield() to file ScalarFieldID_validation/ID_scalarfield.h\n"
]
}
],
"source": [
"# Import the ScalarField.ScalarField_InitialData NRPy module\n",
"import ScalarField.ScalarField_InitialData as sfid\n",
"\n",
"# Output the unity lapse initial data file\n",
"outputname = os.path.join(Ccodesdir,\"outputSFID_unity_lapse-validation.txt\")\n",
"sfid.ScalarField_InitialData(outputname,ID_Family,\n",
" phi0,r0,sigma,NR,RMAX,CoordSystem=CoordSystem,\n",
" sinhA=sinhA,sinhW=sinhW,lapse_condition=\"Unity\")\n",
"\n",
"# Output the \"pre-collapsed\" lapse initial data file\n",
"outputname = os.path.join(Ccodesdir,\"outputSFID_precollapsed_lapse-validation.txt\")\n",
"sfid.ScalarField_InitialData(outputname,ID_Family,\n",
" phi0,r0,sigma,NR,RMAX,CoordSystem=CoordSystem,\n",
" sinhA=sinhA,sinhW=sinhW,lapse_condition=\"Pre-collapsed\")\n",
"\n",
"# Output C codes\n",
"sfid.NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(Ccodesdir=Ccodesdir)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:10.855083Z",
"iopub.status.busy": "2021-06-15T10:05:10.854535Z",
"iopub.status.idle": "2021-06-15T10:05:10.859497Z",
"shell.execute_reply": "2021-06-15T10:05:10.859132Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Unity lapse initial data test: PASSED!\n",
" \"Pre-collapsed\" lapse initial data test: PASSED!\n",
" ADM quantities interpolation file test: PASSED!\n",
" Scalar field interpolation file test: PASSED!\n",
"Scalar field Spherical to Curvilinear test: PASSED!\n",
" Scalar field driver test: PASSED!\n"
]
}
],
"source": [
"import filecmp\n",
"\n",
"if filecmp.cmp(os.path.join(Ccodesdir,'outputSFID_unity_lapse.txt'),\n",
" os.path.join(Ccodesdir,'outputSFID_unity_lapse-validation.txt')) == False:\n",
" print(\"ERROR: Unity lapse initial data test FAILED!\")\n",
" sys.exit(1)\n",
"else:\n",
" print(\" Unity lapse initial data test: PASSED!\")\n",
"\n",
"if filecmp.cmp(os.path.join(Ccodesdir,'outputSFID_precollapsed_lapse.txt'),\n",
" os.path.join(Ccodesdir,'outputSFID_precollapsed_lapse-validation.txt')) == False:\n",
" print(\"ERROR: \\\"Pre-collapsed\\\" lapse initial data test FAILED!\")\n",
" sys.exit(1)\n",
"else:\n",
" print(\" \\\"Pre-collapsed\\\" lapse initial data test: PASSED!\")\n",
"\n",
"if filecmp.cmp(os.path.join(Ccodesdir,'ID_scalarfield_ADM_quantities.h'),\n",
" os.path.join(Ccodesdir,'ID_scalarfield_ADM_quantities-validation.h')) == False:\n",
" print(\"ERROR: ADM quantities interpolation file test FAILED!\")\n",
" sys.exit(1)\n",
"else:\n",
" print(\" ADM quantities interpolation file test: PASSED!\")\n",
"\n",
"if filecmp.cmp(os.path.join(Ccodesdir,'ID_scalarfield_spherical.h'),\n",
" os.path.join(Ccodesdir,'ID_scalarfield_spherical-validation.h')) == False:\n",
" print(\"ERROR: Scalar field interpolation file test FAILED!\")\n",
" sys.exit(1)\n",
"else:\n",
" print(\" Scalar field interpolation file test: PASSED!\")\n",
"\n",
"if filecmp.cmp(os.path.join(Ccodesdir,'ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2.h'),\n",
" os.path.join(Ccodesdir,'ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2-validation.h')) == False:\n",
" print(\"ERROR: Scalar field Spherical to Curvilinear test FAILED!\")\n",
" sys.exit(1)\n",
"else:\n",
" print(\"Scalar field Spherical to Curvilinear test: PASSED!\")\n",
"\n",
"if filecmp.cmp(os.path.join(Ccodesdir,'ID_scalarfield.h'),\n",
" os.path.join(Ccodesdir,'ID_scalarfield-validation.h')) == False:\n",
" print(\"ERROR: Scalar field driver test: FAILED!\")\n",
" sys.exit(1)\n",
"else:\n",
" print(\" Scalar field driver test: PASSED!\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 5: 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.pdf](Tutorial-ADM_Initial_Data-ScalarField.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": 22,
"metadata": {
"execution": {
"iopub.execute_input": "2021-06-15T10:05:10.862225Z",
"iopub.status.busy": "2021-06-15T10:05:10.861827Z",
"iopub.status.idle": "2021-06-15T10:05:13.918364Z",
"shell.execute_reply": "2021-06-15T10:05:13.918792Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created Tutorial-ADM_Initial_Data-ScalarField.tex, and compiled LaTeX file\n",
" to PDF file Tutorial-ADM_Initial_Data-ScalarField.pdf\n"
]
}
],
"source": [
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ADM_Initial_Data-ScalarField\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}