{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Testing Piecewise Polytropic EOS\n", "\n", "## Author: Leo Werneck\n", "\n", "## This notebook assesses the piecewise polytropic Equation of State (EOS) relations used in neutron star models, following [J.C. Read *et al.* (2008)](https://arxiv.org/pdf/0812.2163.pdf) and [Takami *et al.* (2014)](https://arxiv.org/pdf/1412.3240v2.pdf). It evaluates $K$ and $\\Gamma$, formulating a function, and generates relevant plots for testing, concluding with an implementation of the EOS in C code." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Table of Contents\n", "$$\\label{toc}$$\n", "\n", "This module is organized as follows\n", "\n", "1. [Step 1](#introduction): **Introduction**\n", "1. [Step 2](#plugging_in_values): **Plugging in some values**\n", "1. [Step 3](#p_cold): **Taking a look at $P_{\\rm cold}$**\n", " 1. [Step 3.a](#p_cold__computing): *A function to evaluate $P_{\\rm cold}$ based on `IllinoisGRMHD`*\n", " 1. [Step 3.b](#p_cold__plotting): *Plotting $P_{\\rm cold}\\left(\\rho\\right)$*\n", "1. [Step 4](#eps_cold): **Taking a look at $\\epsilon_{\\rm cold}$**\n", " 1. [Step 4.a](#eps_cold__computing): *A function to evaluate $\\epsilon_{\\rm cold}$ based on `IllinoisGRMHD`*\n", " 1. [Step 4.b](#eps_cold__plotting): *Plotting $\\epsilon_{\\rm cold}\\left(\\rho\\right)$*\n", "1. [Step 5](#ppeos__c_code): **The piecewise polytrope EOS C code**\n", " 1. [Step 5.a](#ppeos__c_code__prelim): *Preliminary treatment of the input*\n", " 1. [Step 5.a.i](#ppeos__c_code__prelim__computing_ktab): Determining $\\left\\{K_{1},K_{2},\\ldots,K_{\\rm neos}\\right\\}$\n", " 1. [Step 5.a.ii](#ppeos__c_code__prelim__computing_eps_integ_consts): Determining $\\left\\{C_{0},C_{1},C_{2},\\ldots,C_{\\rm neos}\\right\\}$\n", " 1. [Step 5.b](#ppeos__c_code__eos_struct_setup) *Setting up the `eos_struct`*\n", " 1. [Step 5.c](#ppeos__c_code__find_polytropic_k_and_gamma_index) *The `find_polytropic_K_and_Gamma_index()` function*\n", " 1. [Step 5.d](#ppeos__c_code__compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold): *The new `compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold()` function*\n", "1. [Step n](#latex_pdf_output): **Output this notebook to $\\LaTeX$-formatted PDF file**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 1: Introduction \\[Back to [top](#toc)\\]\n", "$$\\label{introduction}$$\n", "\n", "We will test here a few piecewise polytrope (PP) equation of state (EOS) relations. The main reference we will be following is [J.C. Read *et al.* (2008)](https://arxiv.org/pdf/0812.2163.pdf), but we will also try out [the approach of `Whisky`, presented in Takami *et al.* (2014)](https://arxiv.org/pdf/1412.3240v2.pdf)\n", "\n", "First, we will start with the original one implemented by IllinoisGRMHD:\n", "\n", "$$\n", "\\boxed{\n", "P_{\\rm cold} =\n", "\\left\\{\n", "\\begin{matrix}\n", "K_{0}\\rho^{\\Gamma_{0}} & , & \\rho \\leq \\rho_{0}\\\\\n", "K_{1}\\rho^{\\Gamma_{1}} & , & \\rho_{0} \\leq \\rho \\leq \\rho_{1}\\\\\n", "\\vdots & & \\vdots\\\\\n", "K_{j}\\rho^{\\Gamma_{j}} & , & \\rho_{j-1} \\leq \\rho \\leq \\rho_{j}\\\\\n", "\\vdots & & \\vdots\\\\\n", "K_{N-1}\\rho^{\\Gamma_{N-1}} & , & \\rho_{N-2} \\leq \\rho \\leq \\rho_{N-1}\\\\\n", "K_{N}\\rho^{\\Gamma_{N}} & , & \\rho \\geq \\rho_{N-1}\n", "\\end{matrix}\n", "\\right.\n", "}\\ .\n", "$$\n", "\n", "Notice that we have the following sets of variables:\n", "\n", "$$\n", "\\left\\{\\underbrace{\\rho_{0},\\rho_{1},\\ldots,\\rho_{N-1}}_{N\\ {\\rm values}}\\right\\}\\ ;\\\n", "\\left\\{\\underbrace{K_{0},K_{1},\\ldots,K_{N}}_{N+1\\ {\\rm values}}\\right\\}\\ ;\\\n", "\\left\\{\\underbrace{\\Gamma_{0},\\Gamma_{1},\\ldots,\\Gamma_{N}}_{N+1\\ {\\rm values}}\\right\\}\\ .\n", "$$\n", "\n", "Also, notice that $K_{0}$ and the entire sets $\\left\\{\\rho_{0},\\rho_{1},\\ldots,\\rho_{N-1}\\right\\}$ and $\\left\\{\\Gamma_{0},\\Gamma_{1},\\ldots,\\Gamma_{N}\\right\\}$ must be specified by the user. The values of $\\left\\{K_{1},\\ldots,K_{N}\\right\\}$, on the other hand, are determined by imposing that $P_{\\rm cold}$ be continuous, i.e.\n", "\n", "$$\n", "P_{\\rm cold}\\left(\\rho_{0}\\right) = K_{0}\\rho_{0}^{\\Gamma_{0}} = K_{1}\\rho_{0}^{\\Gamma_{1}} \\implies\n", "\\boxed{K_{1} = K_{0}\\rho_{0}^{\\Gamma_{0}-\\Gamma_{1}}}\\ .\n", "$$\n", "\n", "Analogously,\n", "\n", "$$\n", "\\boxed{K_{j} = K_{j-1}\\rho_{j-1}^{\\Gamma_{j-1}-\\Gamma_{j}}\\ ,\\ j\\in\\left[1,N\\right]}\\ .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 2: Plugging in some values \\[Back to [top](#toc)\\]\n", "$$\\label{plugging_in_values}$$\n", "\n", "Just so that we work with realistic values (i.e. values actually used by researchers), we will implement a simple check using the values from [Table II in J.C. Read *et al.* (2008)](https://arxiv.org/pdf/0812.2163.pdf):\n", "\n", "| $\\rho_{i}$ | $\\Gamma_{i}$ | $K_{\\rm expected}$ |\n", "|------------|--------------|--------------------|\n", "|2.44034e+07 | 1.58425 | 6.80110e-09 |\n", "|3.78358e+11 | 1.28733 | 1.06186e-06 |\n", "|2.62780e+12 | 0.62223 | 5.32697e+01 |\n", "| $-$ | 1.35692 | 3.99874e-08 |" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:05.097985Z", "iopub.status.busy": "2021-05-17T23:43:05.097085Z", "iopub.status.idle": "2021-05-17T23:43:05.246891Z", "shell.execute_reply": "2021-05-17T23:43:05.247383Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "K_ppoly_tab[1] = 1.0618444833278535e-06 | K_expected[1] = 1.06186e-06 | Diff = 1.5516672146434653e-11\n", "K_ppoly_tab[2] = 53.2750845839615 | K_expected[2] = 53.2697 | Diff = 0.005384583961500766\n", "K_ppoly_tab[3] = 3.999206917243191e-08 | K_expected[3] = 3.99874e-08 | Diff = 4.669172431909315e-12\n" ] } ], "source": [ "# Determining K_{i} i != 0\n", "# We start by setting up all the values of rho_{i}, Gamma_{i}, and K_{0}\n", "import numpy as np\n", "rho_ppoly_tab = [2.44034e+07,3.78358e+11,2.62780e+12]\n", "Gamma_ppoly_tab = [1.58425,1.28733,0.62223,1.35692]\n", "K_ppoly_tab = [6.80110e-09,0,0,0]\n", "K_expected = [6.80110e-09,1.06186e-06,5.32697e+01,3.99874e-08]\n", "NEOS = len(rho_ppoly_tab)+1\n", "\n", "for j in range(1,NEOS):\n", " K_ppoly_tab[j] = K_ppoly_tab[j-1] * rho_ppoly_tab[j-1]**(Gamma_ppoly_tab[j-1] - Gamma_ppoly_tab[j])\n", " print(\"K_ppoly_tab[\"+str(j)+\"] = \"+str(K_ppoly_tab[j])+\\\n", " \" | K_expected[\"+str(j)+\"] = \"+str(K_expected[j])+\\\n", " \" | Diff = \"+str(np.fabs(K_ppoly_tab[j] - K_expected[j])))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 3: Taking a look at $P_{\\rm cold}$ \\[Back to [top](#toc)\\]\n", "$$\\label{p_cold}$$\n", "\n", "\n", "\n", "## Step 3.a: A function to evaluate $P_{\\rm cold}$ based on `IllinoisGRMHD` \\[Back to [top](#toc)\\]\n", "$$\\label{p_cold__computing}$$\n", "\n", "The results above look reasonable to the expected ones. Let us then see what the plot of $P_{\\rm cold}\\times\\rho$ looks like. For the case of our specific example, which has $N_{\\rm EOS}=4$ (where $N_{\\rm EOS}$ stands for the number of polytropic EOSs used), we have\n", "\n", "$$\n", "\\boxed{\n", "P_{\\rm cold} =\n", "\\left\\{\n", "\\begin{matrix}\n", "K_{0}\\rho^{\\Gamma_{0}} & , & \\rho \\leq \\rho_{0}\\\\\n", "K_{1}\\rho^{\\Gamma_{1}} & , & \\rho_{0} \\leq \\rho \\leq \\rho_{1}\\\\\n", "K_{2}\\rho^{\\Gamma_{2}} & , & \\rho_{1} \\leq \\rho \\leq \\rho_{2}\\\\\n", "K_{3}\\rho^{\\Gamma_{3}} & , & \\rho \\geq \\rho_{2}\n", "\\end{matrix}\n", "\\right.\n", "}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:05.252941Z", "iopub.status.busy": "2021-05-17T23:43:05.252228Z", "iopub.status.idle": "2021-05-17T23:43:05.254115Z", "shell.execute_reply": "2021-05-17T23:43:05.254589Z" } }, "outputs": [], "source": [ "# Set the pressure through the polytropic EOS: P_cold = K * rho^{Gamma}\n", "def P_cold(rho,rho_ppoly_tab,Gamma_ppoly_tab,K_ppoly_tab,NEOS):\n", " if rho < rho_ppoly_tab[0] or NEOS==1:\n", " return K_ppoly_tab[0] * rho**Gamma_ppoly_tab[0]\n", "\n", " if NEOS > 2:\n", " for j in range(1,NEOS-1):\n", " if rho_ppoly_tab[j-1] <= rho and rho < rho_ppoly_tab[j]:\n", " return K_ppoly_tab[j] * rho**Gamma_ppoly_tab[j]\n", "\n", " if rho >= rho_ppoly_tab[NEOS-2]:\n", " return K_ppoly_tab[NEOS-1] * rho**Gamma_ppoly_tab[NEOS-1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.b: Plotting $P_{\\rm cold}\\left(\\rho\\right)$ \\[Back to [top](#toc)\\]\n", "$$\\label{p_cold__plotting}$$\n", "\n", "To create a reasonably interesting plot, which will also test all regions of our piecewise polytrope EOS, let us plot $P_{\\rm cold}\\left(\\rho\\right)$ with $\\rho\\in\\left[\\frac{\\rho_{0}}{10},10\\rho_{2}\\right]$.\n", "\n", "**Remember**: We *don't* expect $P_{\\rm cold}$ to be *smooth*, but we *do* expect it to be *continuous* (by construction)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:05.264933Z", "iopub.status.busy": "2021-05-17T23:43:05.264135Z", "iopub.status.idle": "2021-05-17T23:43:06.146120Z", "shell.execute_reply": "2021-05-17T23:43:06.146722Z" } }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAEiCAYAAAAoKdyhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABPVUlEQVR4nO3dd3wUdf7H8dcHQu810kFAinRCKDasZzuxYEORKuCdXU/9nfU89Q71bKd4Ir2DeKeeBRQ1ViAhECD03gk9JED65/fHDOcaQ7JJNju7yef5eOwju7OzM++dTb757Mx3viOqijHGGGOMKXnlvA5gjDHGGFNWWOFljDHGGBMkVngZY4wxxgSJFV7GGGOMMUFihZcxxhhjTJBY4WWMMcYYEyRWeBljjDHGBIkVXsYYY4wxQWKFVwgRkaEioj63FBFZKSL3ikiEO89zIlLoUW9F5HoRebgY2QaLyE6fx2tF5A+55mkqIv8UkcUictJ9Dy2Luk6f5Y7MtV3SRCRRRIYUd9nGeKkU/M0PFJEPRWSHiJwSkQ0i8jcRqVGM9Rbp/QaTiLwlIp/mmuZ5OyUiD4rIahEp0f/t4foZBWi5xd7GVniFppuBvsBNQCzwT+CZYi7zeqDIjTDQE4gHEJHqQLvTj320AW4BjgI/FGNduXUH0nC2SV/gBuA4MEVELg7geozxSrj+zT8KZAN/Bq4E3gXuAb4q6X/+XhGR1sAY4LlcT4VCO/Ue0AAo019K8/mMAqHY2zgicFlMACWo6mb3/pci0gZ4gOI3xMXRE1jo3u8B5AArc83zvapGgvPtD7giQOvuBqxV1SWnJ4jIdmAtcDXwbYDW4xcRaQS0UtWf83iuKnCRqn4RzEwm7IXr3/zvVfWgz+PvROQIMBXoD3xT0iE98CCwUlWX5ZreDY/bKVU9JSLTcAriyWearwy0YQ+S92dUbP5u4/yUym8kpVAcUFNEGp5pBhG50j3Ed0pEkkXkIxFp5z43Bac6b+KzG3y7vyt3v7l245dvu1E4DUya73yqmlOYN+XnugXoAqzO9dRx92fVQK/TDw8CC0XkfN+JIlIF+C8w2d1DkC8RiRCRR0Vkq4hkiMhBEXmuRBKbcBMuf/MH+a0492cTf9fnR54zvtdc890uIuvdw3yrReQ6EYkRkZgA5agE3AnMyjU9lNqpOUBHEemXzzwPEoA2LNdrQ/oz8nk+EO2uP9v4jKzwCg+tcHbnp+b1pIhcCXzmPn8rzq7+TsCPItIE+CvwOXCQX+8Gz5eIbHeP42cD1YHP3cf/ALr4NOgtC/NmRKSl+7rn/Ji9rbvuVbmmX+T+DPg3Gj88CSwAvhCR8+BXDVZX4ApVzfOzOk1EKrjzPwj8HbgOiAGeFZFbSyy5CRfh/Dd/+m9znc9yC/M3nztTQe/19HyXAzOB9cCNwKvAG8A5hV1nPvoAtfltV4pQaqcSgBScQ79nUuw2zFeYfEaBbHcTKHgbn5mq2i1EbsBQQHH6UkQAdYDROI3gR+48zzkf269etwzYBET4TGsFZAKvuY+nALsLmacjzrfe14A17v1uON/iHvJ5XDGP145030vLPJ5rAWQBz/iR4RZ3OVf6bJMbgSSchr2yR59VBDDf3RaXAouAQ0BXP1//BJAMNPeZVsFd3iyvfxftFpxbafqbd1/fBDgAfJVremH+5n/1fv15r+60n4FEQHym9XS3b0yAPq/HcQ65Vsw1PaTaKZyi48sC5ilyGxaOn5H7XMDaXX+28ZlutscrNK3H+YU9AozD+YYwPK8ZRaQaTv+LuaqadXq6qm4DfuKXb1yFpqprVTUBaIbzR5EAnABqAB+oaoJ7yyjkcneoaoSqPu/H7N3cn1/wyzaZg/Mt5WLNdejDX+6379r5PH8sv2/17ra+DaexWoTzGVymqrn7wOS17HI4nZ7fUdX/nTWmqpnADpxvzohIVTdHeb/elAlnYf837x6a+hinwBqWa7mF+Zv3XaZf79X9G4kCPlT3v6I7XzywLdcynxaRjSKSIyLX53qutYj86D6/QkSickVqDBzP4/13c38GtJ1yM52xrcqnnTroZj2j4rRhuTKExWdUAu1ugdv4TMKm8BKRSSJyQEQS/Zj3QhFZLiJZIjIw13NDRGSTewvVMz9uAHoB7YFqqnqXqh45w7x1AAH25fHcfqBuUQKISHn3WHgEcB6w2L1/AbAH2O8+L0VZfiF0B3bjbI8o4Fyglqreqqr7S3jdBYnA+WPNwfnW5G8/jmics2I+853obssm/NII9cDpdJ0dkLQmlIX137zPYaqzgd+p6u6iZMiDv++1Ps7f4IE85kvK9fgrnD1T3+cx77+Aqap6DvAYMDPX+60MpOfxulBrp04BVfyYr6htmK9w+YwC3e76u41/I2wKL5zd5v4eT92Jsws/dwfIusCzQG+cD+FZEakTuIgBk6iqy1R1gx/flI7i7KY9K4/nzsL55lUUX+N8c8sEGgHT3fsTcX5JTz9X5G/XfuoGLHO3R7z7jfyU7wwi0ltEvheRVe7tDp/nLheRZe70n0SkR14rEZEB4nT4XCkiYwsKJSKVgY9w/kj74DQUC0Skrx/vqY/7M/c/p344jdhH7uNoYK+I/Fec8ZG+DtHfV1N8Yfs37/abmY9TcFytqrk7mBeHv+/1kJstr5MRIn0fqOoSVd2aeyYRaYDztznFne8rnIKip89sh3H6D+XWjQLaKXcdebZVJdBO1cXZJmdUzDbMV7h8RoFudwvcxmcSNoWXqn5PrgbF3eW4QETiReQHEWnvzrtdVVfhVPG+fofT9+CIqh7ll6o6bKnqCZwzj2723TUqIi1wfqFi3EnpFK46H43z7e1VYLN7vxfO7tWnfB7nHtcnYEQkEucPd0U+89TBObzxjKp2wekYutB9rgHO7v6R7nMvAf8RkYq5ltEQ57TggaraFdgK1MpnnZWA/+C8/8tUNQ6nQ+kinIarz5le6zrdSLT2WWY54G84nXO/cyf3ApoCg3D2hBwDfjWApSl7Qulv3v29nQlcAlyvPkMpBIK/79XdO7EMuMl374eI9MTpa+SP5sA+99DTadvd6aetByqKSFOfdRTYTrnz5dlWlVA71QrYkE+W4rZh/xMOn5Er0O1uvts4P2FTeJ3BeOA+Ve2JM6bGuALmbwLs8nm8mwCe8uyhp3HOqvlURH4vIrfjFJXJOGcjgTOWTF0RuUdEeolI5/wW6H7zXoZztsln7v0UnN3FE91vdstUNcX3deKMZD2QX37Jr3KnXeQzTwv3MHBBYxR1d3/m16D1AzapaoybW1X19LeQPsBqt58Kqnp6F3PuM2j6AKtU9fRh7Ik4nZvP5EWcb0WX+Sw7E6fh+hr4r+R/KnZPnN+9cSJyk4jc6L6uO3C7/jIsRzRwv6qmuH0i4sj1zdCUWaHyN/8OzuCv/wBOiEgfn5tvceLv33xR3ys4RzPOxSlarhaRu4APcA53BWqom9OHvqJ9pvnTTsGZ26qAtlPi9Ac7h7wP051W3DYst1D/jCCA7a6f2/iMwrbwcn8p+gEfiEgCzmiyjTwN5RFVXQBcg7N7dR7OMfB1wPmqutedbQLOt6qXcEbG/m9By3W/cV2Kc9oxwFXAigL6LHzg3sa4j8e5j//iu2igPAX//nVzfxbUoBWGv5e5yG++F4H+qvqrXD4N13V6hlOxxemI2s5dxkc4v7fTcfoL9FXVte589XB+n307uvbDm+EzTIgJob/5q9yfTwKLc91G+i4a//7mf8PP93r6sNMdQAecvTmPA4/g/FNP9mNVO4FG7qHT01q600+vYzvOtvy9zzzd3J+h0k5dA2TgbIMzKXIblmeIEP+MSqDd9Wcbn5kG8RTX4t7cDZzo3q+Js8sxv/mn4OyWPf34duA9n8fv4VS6nr83uxX5d6IOzh9tf/exAPXd+w1wjsF3cR9fhbNbuiJOg1XbZ77DwLnu47s5w1AYAch7vrvsngXMdyXON8D27uMbcAZnzPM0frvZzW6/veEcMkoDns7juRicw6O+074G7nbvX44zRILkmmcoTpFQtZBZ8myr8mun3MeKU9D41U7hnFk53ettH0qfUaDb3eJuY883eiE/oJa4hZf7+GfgZve+kGv8EX5beNXFOXOhjnvbBtT1+n3Zrdi/F72BH3GO068CBvk8dzlO/4OVOKc293Cn/6/wch8PwOkbsBJnYL1juRu0AGV9AKfvTb4FFM6u+/dwLjOyBvgUaOL1trab3UL1htOf7V2c611ehDOkxTqcvsGNfOZ7DueQU7pb8OwGmrrPtXX/r2zEGSQzOo/1RLjLfbQIGfNsq87UTrnP+X5JzLedwtn7lg608frzCKXPKJDtbiC2sbgLCnkiMhvn2l/1cU49fRbnOmDv4uwarADMUdXnRaQXzi7AOjiV9H5VPdddznCcC7oCvKiqk4P5PkzZJs41vs5Vp1+iMSZA3MOkc3H6QtXDGX/sB+DP+ku/qECtqw9OcVRQv+KgEmcE+TqqOtvrLHnx6jMKZLsbiG0cNoWXMcYYY0y4C9vO9cYYY4wx4SbC6wD+qF+/vrZs2dLv+U+cOEG1atVKLlAAhVNWsLwlzfL+Ij4+/pCqNiiRhQdRYduvQAuF3ynL4P36LUNwM+Tbfnnd2c6fW8+ePbUwvv3220LN76VwyqpqeUua5f0Fzkjgnrc/xb0Vtv0KtFD4nbIM3q/fMgQ3Q37tlx1qNMYYY4wJEiu8jDHGGGOCxAovY4wxxpggscLLGGOMMSZIrPAyxhhjjAkSK7yMMcYYY4LECi9jjDHGmCCxwssYU+K+3XCAObE7vY5hjDGFdvREBs99soa0zOyALM8KL2NMiVq79zj3zlzO9CU7yMjK8TqOMcb4bV/yKW5+bzGzYneyZm9yQJYZFpcMMsaEp/3JaQyfEkeNyhWYOKQXFSPsu54xJjxsPZjK4ImxJJ/KZNrwaHq2qBuQ5VrhZYwpESlpmQybEkdKWiYfjOnHWbUqex3JGGP8krgnmSGTYgGYM6oPnZrUCtiyrfAyxgRcZnYOf5y1go1JKUwa2ouOjWt6HckYY/yyeMth7p62jFpVKjB9RDRnN6ge0OVb4WWMCShV5emPEvl+40H+fmNnLjqngdeRjDHGL1+u2c+9s1fQom5Vpo2IplGtKgFfhxVexpiAGhezhTlxu/jjxa25Lbq513GMMcYv8+N38/iHq+jUpBZThvaiTrWKJbIeK7yMMQHzccIeXlm4gQHdGvPoFe28jmOMMX6Z8MNWXvhsHee3qc97g3tSrVLJlUdWeBljAmLJ1sP86YNV9G5Vl5cHdkFEvI5kjDH5UlVeWbiBcTFbuLrzWbx+azcqRZQv0XVa4WWMKbbNB1IYNW0ZzepWYfzgqBJvuIwxpriyc5SnPkpkduxObo9uzgvXd6J8uZL/wmiFlzGmWA6mpDN0chwVI8oxZVg0tapW8DqSMcbkKz0rm4fnruSz1fv4Q//W/Ol37YK2l94KL2NMkaVnKSOmxnE4NYM5o/rQrG5VryMZY0y+0rKUkVOX8cOmQzx5dQfuvvDsoK7fCi9jTJFk5yjvrkwn8dBJ3hscRddmtb2OZIwx+Tp6IoOX49LYkXKKVwZ24eaoZkHP4Mn1O0SktojMF5H1IrJORPp6kcMYUzSqyvP/XUPCwWye/f25XN4x0utIxhiTr/3Jadzy3mJ2puTw7h09PCm6wLs9Xm8CC1R1oIhUBOz4hDFhZOKP25i6eAe/axnBkH4tvY5jjDH58r3u4iM9K3PFuWd5liXoe7xEpBZwITARQFUzVPVYsHMYY4rmi9X7ePHzdVzV6SxubVcyAwwaY0ygJO5J5uZ/LSYtM5s5o/rQoZ63Z117scerFXAQmCwiXYF44AFVPeE7k4iMAkYBREZGEhMT4/cKUlNTCzW/l8IpK1jekhbqeTcfzWZsXBqta5XjhkbHOXniREjnLQ4RmQRcCxxQ1U7utLrAXKAlsB24RVWPepXRGJO/JVsPc/fUZdT0ue5izCZvM3lReEUAPYD7VHWpiLwJPAE87TuTqo4HxgNERUVp//79/V5BTEwMhZnfS+GUFSxvSQvlvNsPneDhd3+mcZ2qzL2nH/WqVwrpvAEwBXgbmOYz7Qnga1X9u4g84T5+3INsxpgCfLU2iT/OWk7zulWZXkLXXSwKLzrX7wZ2q+pS9/F8nELMGBOiDqemM3RyLKrKlGHR1KteyetIJU5VvweO5Jo8AJjq3p8KXB/MTMYY/3wYv5sxM+LpcFYN5o3uGzJFF3iwx0tV94vILhFpp6obgEuBtcHOYYzxz8mMLIZPXca+5DRm3d2HVvWreR3JS5Gqus+9vx/I83TO4nSVCLRQOHxtGbxff1nKsHB7JrPXZ3BuvXKMaZ/Jqrifg54hP16d1XgfMNM9o3ErMMyjHMaYfGRl53DvrBWs3n2Mf93Zk54t6ngdKWSoqoqInuG5IneVCLRQOBxsGbxff1nIoKr848uNzF6/Od/rLnq9HTwpvFQ1AYjyYt3GGP+oOtcx+2b9AV68oZOnp1+HkCQRaaSq+0SkEXDA60DGGGdA56c/TmTW0p3cHt2MF67vHJTrLhaFJwOoGmNC35tfb2JO3C7uu6QNd/Ru4XWcUPEJMMS9PwT42MMsxhggIyuH++esYNbSnfyhf2teuiF0iy6wSwYZY/IwJ3YnbyzaxMCeTXn48nO8juMJEZkN9Afqi8hu4Fng78A8ERkB7ABu8S6hMeZEehZjZsR7dt3ForDCyxjzK1+vS+LJjxK56JwG/O3GzoiE7jfHkqSqt5/hqUuDGsQYk6ejJzIYNiWOVbuP8fLALtzi0SWACssKL2PM/6zYeZQ/zlpOx0Y1GXdHDyqUt94IxpjQsz85jcETl7LjyEnevbMnvwujPqhWeBljANh26AQjpi6jYY3KTBrai2qVrHkwxoSebYdOcOeEpSSfymTKsF70a13f60iFYi2rMYaDKekMmRQLwNTh0TSoUfoHSDXGhJ/EPckMnRxLjsLsu/vQuWktryMVmhVexpRxJ9KzGDE1jgMpacy2AVKNMSFq6dbDjHSvuzhtRDStG1T3OlKRWOFlTBmWmZ3DH2ctJ3FPMu/fFUX35jZAqjEm9Cxyr7vYrG5Vpg2PpnHt0LkEUGFZ4WVMGaWq/Pnfq4nZcJC/3diZSzvkefUbY4zxTHaOMjduF09/nEinxjWZPCyautUqeh2rWKzwMqaMev2rjXwQv5sHLm3L7dHNvY5jjDH/c+xkBnPjdjF9yQ52Hz3F+W3q86/BPaleCk76Cf93YIwptJlLd/DWN5u5NaoZD17W1us4xhgDwLp9x5n683Y+SthDWmYOvVvV5cmrO3B5x0giSsnwNlZ4GVPGLEjcx9MfJXJxuwa8eEOnMjtAqjEmNGRl5/DV2iSm/LydpduOULlCOW7o3oS7+rakQ6OaXscLOCu8jClDFm85zP1zEujarDbv3NGj1HyDNMaEnyMnMpgdu5OZS3awNzmNpnWq8Oer23NLVDNqVw3vflz5scLLmDJizd5kRk1bRou6VZk8tBdVK9qfvzEm+BL3JDPl5+18snIvGVk5nNemHs9ddy6XdogM6YtbB4q1vMaUATsOn2DIpDhqVI5g2ojoUv1t0hgTejKzc1iQuJ83l5xi84IfqVqxPLdENWVI35a0jazhdbygssLLmFLuYEo6d02KJSsnhzmj+tKoVviOf2OMCS8HU9Kdw4lLd5B0PJ2GVYWnr+3IwJ5NqVWlgtfxPGGFlzGl2PG0TIZMiuXA8XRm3d2bNg3L1jdLY4w3EnYdY+rP2/ls1T4ysnO48JwG/O3GFrBvLZec38rreJ6ywsuYUiotM5tR05axMSmFCUNsVHpjTMlKz8rm89X7mPLzDlbuOkb1ShEM6t2cwX1b/O/yPjH713mc0ntWeBlTCmXnKA/OSWDJ1iO8cWs3+rdr6HUkY0wplXQ8jZlLdjArdheHUtM5u0E1/nLdudzUs2mpGPA00GyLGFPKqCpPf5zIgjX7efrajlzfvYnXkYwxpYyqsnznUab8vIMvVu8jW5VL2jVkSL+WnN+mPuXKwNmJRWWFlzGlzOuLNjFr6U7u6d+aEWW8L4UxJvDith/hhU/XsnJ3MjUqRzCkX0vu6tuCFvWqeR0tLHhSeInIdiAFyAayVDXKixzGlDZTf97OW19v4paopjz2u3ZexzHGlDKTftzGXz9bS6OalXnh+k7c0L0J1exwYqF4ubUuVtVDHq7fmFLl01V7ee6/a7isQyQv3dDZLgVkjAmYnBzlhc/WMemnbVzRMZI3butmgzAXkW01Y0qBHzcd4qG5CUS1qMPbg7rbpYCMMQGTlpnNQ3MT+CJxP0P7teTpazuWiRHmS4pXhZcCX4qIAu+p6vjcM4jIKGAUQGRkJDExMX4vPDU1tVDzeymcsoLlLWlFybstOZuxsWlEVhGGtk5nyU8/lEy4PITb9jXGFM6RExncPW0Zy3ce5alrOjDi/Fa2N72YvCq8zlfVPSLSEPhKRNar6ve+M7jF2HiAqKgo7d+/v98Lj4mJoTDzeymcsoLlLWmFzbv5QAoP/Wsx9WtW4cN7+hFZs3LJhctDuG1fY4z/dhw+wdDJcew5dop3BvXg6s6NvI5UKnhSeKnqHvfnARH5DxANfJ//q4wxvnYdOcmdE2IpX64cM0b0DnrRZYwpvVbsPMrIqcvIVmXWyN5EtazrdaRSI+gdQUSkmojUOH0fuAJIDHYOY8LZgZQ0Bk9cysmMLKaPiKZlfTuN2xgTGAvX7Of295dQrVIE/76nnxVdAebFHq9I4D/uMeIIYJaqLvAghzFhKflUJndNjCXpeDozRvamQ6OaXkcyxpQSU37axl8+XUuXprWZOCSK+tUreR2p1Al64aWqW4GuwV6vMaXByYwshk+JY8vBVCYO6UXPFnb9RWNM8eXkKH/7Yh3v/7CNyztG8tZt3alSsbzXsUolG07CmDCRkZXDmBnLWbHzKG8P6sGF5zTwOpIxphRIy8zmkXkr+Wz1Pu7q24Jnf3+uDRdRgqzwMiYMZOcoD81N4PuNBxl7U2c7u8gYExBH3eEilu04ypNXd2DkBTZcREmzwsuYEKeqPPmf1Xy2eh9PXt2BW3s19zqSMaYU2Hn4JEMnx7L72CneHtSda7s09jpSmWCFlzEhTFX5+xfrmRO3i3svbsPdF57tdSRjTCmQsOsYI6bEka3KzJG96WVnLgaNFV7GhLBxMVt47/utDO7TgkeuOMfrOMaYUuCrtUncN3s5DWpUYsqwaFo3qO51pDLFCi9jQtT0JTt4ZeEGBnRrzF+uO9f6XRhjim364u08+8kaOjWpxcQhvWhQw4aLCDYrvIwJQR8n7OGZjxO5pH1DXr25K+XsDCNjTDHk5ChzN2TwxbY1XNahIW/d3p2qFa0E8IJtdWNCzDfrk3hk3kp6tazLuDt6UKF80C8wYYwpRdIys3n0g5V8sS2TwX1a8Nx1NlyEl6zwMiaErDuczRuLltO+UQ0mDomicgUbwNAYU3THTmYwalo8sduPcEu7Cjw/wLoteM0KL2NCRPyOo7yxPI3m9aozdVg0NSpX8DqSOQMReQgYCSiwGhimqmnepjLm13YdOcmQybHsPnKKt27vTs2jG63oCgF2DMOYEJC4J5mhk2OpVUmYObI39ez6aCFLRJoA9wNRqtoJKA/c5m0qY35t1e5j3DDuZw6lpDN9RDTXdbUxukKF7fEyxmObklK4a1IsNSpF8HC3CjSsWdnrSKZgEUAVEckEqgJ7Pc5jzP98vS6Je2etoF71iswZ1Zs2DWt4Hcn4sMLLGA9tP3SCOyYspXw5YdbdfdieGOd1JFMAVd0jIq8CO4FTwJeq+qXvPCIyChgFEBkZSUxMTNBznpaamurp+i1DcNf/zc5Mpq/NoEXNcjzYDXavjWf32uBmyI9lsMLLGM/sPnqSOyYsJTM7h7mj+9KyfjW2ex3KFEhE6gADgFbAMeADEblTVWecnkdVxwPjAaKiorR///4eJHXExMTg5fotQ3DWn5OjvLxwA9PWbuGS9g355+3dqVbp1//ivd4GlsFhhZcxHjhwPI07JyzleFoms+/uwzmRdiggUESkEtAYqAIcVNWDAV7FZcC208sVkX8D/YAZ+b7KmBKSnpXNnz5YxScr9zKod3Oev+5cImwYmpBlhZcxQXbkRAZ3TFjKgZR0po/oTacmtbyOFPZEpAZwJ3A7EA1UAARQEdkDLATGq2ogjuXuBPqISFWcQ42XAssCsFxjCi35ZCajpi9j6bYjPH5le8ZcdLaduRjirPAyJoiST2UyeOJSdh45yeRhvejZoo7XkcKeiDwMPAlsBT4BXsTp7H4KqAt0Ai4AvhKRJcB9qrqpqOtT1aUiMh9YDmQBK3APKxoTTLuOnGTYlDh2Hj7Jm7d1Y0C3Jl5HMn6wwsuYIDmRnsWwybFsTEph/OAo+rWu73Wk0qIPcJGqJp7h+VhgkojcAwwHLgKKXHgBqOqzwLPFWYYxxbF6dzLDp8aRnpnNtBHR9Dm7nteRjJ+s8DImCNIysxk5dRkJu47xzqAeXNy+odeRSg1VvcXP+dKAcSUcx5gS9+36A/xx1nLqVK3IrJG9aWt9RMOKFV7GlLCMrBzumRHPkm2Hee2WrlzVuZHXkUo9Ebk6v+dV9fNgZTEmkGYt3cnTHyfSoVENJg3pZeP+hSErvIwpQVnZOTwwZwXfbjjISzd05obuTb2OVFbc7P6MBPoCX+N0tr8YWAxY4WXCiqry6pcbeOfbLVzcrgFvD+rxm+EiTHiwT82YEpKVncND81byReJ+nr62I4N6N/c6UpmhqsMARGQR0EFV97uPz8KGfTBhJiMrh8fmr+SjhL3cHt2Mvw7oZMNFhDHPCi8RKY9zCvYeVb3WqxzGlITsHOVP81fx35V7eeKq9ow4v5XXkcqqpsAhn8eH3WnGhIXkU5mMmR7P4q2H+dPv2vGH/q1tuIgw5+UerweAdUBNDzMYE3A5OcrjH67iPyv28OgV5zDmotZeRyrL5gA/ich/AAWuB2Z7msgYP+05doqhk2LZfvgEb9zajeu723ARpYEn+ypFpClwDTDBi/UbU1JycpQnP1rN/PjdPHBpW+69pK3Xkco0VX0OuBdnTK804F5V/YunoYzxw6akFG4a9zP7j6cxdXi0FV2liFd7vN4AHgPOeA5scS4y6/UFMAsjnLKC5c2PqjJ9bQbf7Mri2rMr0C1iDzExewu1DNu+geeOVm9XHzdhI2HXMYZOjqVC+XLMG92XDo3swFBpEvTCS0SuBQ6oaryI9D/TfMW5yKzXF8AsjHDKCpb3TFSVv/x3Ld/s2s7oC8/miavaF6kfhm3fwBCROJxDi795ClBVjQ5yJGP88tPmQ9w9bRn1q1dixojeNK9X1etIJsC82ON1HnCdO85OZaCmiMxQ1Ts9yGJMsakqL32+jik/b2f4ea2KXHSZgBrodQBjCmtB4j7un51Aq/rVmD4i2sboKqWC3sdLVf9PVZuqakvgNuAbK7pMuFJVXl64gfd/2MaQvi14+toOVnSFAFXdcfoGpANd3FuaO82YkDI3bid/mLmcTk1qMm90Xyu6SjEbCMSYYnj9q428G7OFQb2b89x151rRFWJEZBDwI87JPNcCP4jIbd6mMubX3vtuC49/uJrz2zZgxsje1KpawetIpgR5OoCqqsYAMV5mMKao3vp6E299s5lbo5rxwoBOVnSFpseBXqp6FEBE6uC0OXO8DGUM/LLH/N2YLVzbpRGv3dKNihG2P6S0s5HrjSmCcTGbee2rjdzUoyl/u7Ez5cpZ0RWiygGpPo9TsT39JgRk5yhPfbSa2bG7uKN3c54f0Iny1o6UCVZ4GVNI47/fwssLNnB9t8a8PLCLFV2hbQbws4h86D6+EZjuYR5jSM/K5uG5K/ls9T7uvbgNj1xxju0xL0Os8DKmEN77bgt/+2I913RpxKs3d7VvqCFOVceKyNc4Z1MD3KOq8V5mMmXbifQsxsyI54dNh3jqmg6MvOBsryOZILPCyxg/vRuzhbEL1vP7ro15/ZaudpHaMCAiLwMvqeoy93EdEfm7qj7hcTRTBh07mcHQyXGs2n2MVwZ24eaoZl5HMh6w/xzG+OGdbzczdsF6rrOiK9xcrqrHTj9wO9lf4V0cU1YdTcvhlvcWs3bfcd69s6cVXWVYofd4uQOfnpGqfl70OMaEnre/2cSrX27k+m6NefVmK7rCTDkRqaGqKQAiUhOwc/VNUG0/dIIXl6ZxKrscU4b1ol/r+l5HMh4qyqHGm92fkUBf4Gucy3BcDCwGrPAypcZbX2/ita82cmP3JrxifbrC0ZvAjyIy1318K/C6h3lMGbN273HumhRLepYye3QfujSt7XUk47FCF16qOgxARBYBHVR1v/v4LJwziIwpFd5YtJE3Fm3ixh5NeGWgFV3hSFUniUgszhdDgEGqusbLTKbsiNt+hOFT4qheKYKHe1exossAxetc3xQ45PP4sDvNmLCmqry+aBNvfb2JgT2bMvamLlZ0hQER6Qr0AJKARaqaAaCqiUCil9lM2fPt+gPcMzOexrWrMH1EbzYlLPU6kgkRxSm85gA/ich/AAWuB2YHIpQxXlFVXv9qI299s5lbopry9xttnK5wICKjgHdxuj0AbBKRS1R1j4exTBn1ccIeHpm3kvaNajB1WDT1qldik9ehTMgoci9hVX0OuBc4BaQB96rqXwKUy5igU1X+8aVTdN3Wq5kVXeHlMWAccBbQCzgAjPU0kSmTpi3ezoNzE+jZog6z7+5DveqVvI5kQkyxxvFS1TggLkBZjPGMqvLKwg2Mi9nC7dHNePF6uwxQmGkBvKqqB4ADIjIUWO1tJFOWqCpvfb2Z1xdt5LIOkbw9qDuVK5T3OpYJQUUZTiIO59Dib54CVFWji53KmCBSVcYu2MC/vtvCoN7NeWFAJyu6wk95nL3vAKjqFhFBRBqp6j4Pc5kyICdHef7TtUz5eTs39WjK2Js627Az5oyKssdrYMBTGOMRVeWvn65j0k/buLNPc56/zoquMDZKRH4GElT1CJANVPE4kynlMrNzeHz+Kv69Yg/Dz2vFU9d0sDbE5Ksow0nsOH3fHUKil/swVlWTAhXMmJKWk6M89XEis5buZNh5LXnm2o52odrw9S3wMPA8oCKyF6iMU4x9DSxzR603JmDSMrO5d9ZyFq07wKNXnMMfL25jbYgpUJH3hYrIIOBH4BrgWuAHEbktUMGMKUnZOcpjH65i1tKd3NO/tRVdYU5VL1XVukAb4DZgJk4xNhJYCBwSETuxzATM8bRM7poUy9frD/DX6ztx7yVtrQ0xfilO5/rHgV6nv0WKSB0gBmeYCWNCVmZ2Dg/PW8l/V+7locvO4f5L7VtqaaGqW4GtwAenp4lISyAKZ4wvY4rtUGo6QybFsmF/Cm/e1p3rujb2OpIJI8UpvMoBqT6PU7GLbpsQl56Vzf2zV7BwTRJPXNWeMRe19jqSKSYRaaWq2870vKpuB7YD88WpsJuq6q4gxTOlzO6jJ7lrYix7k0/x/pAoLm7X0OtIJswUp1CaAfwsIk+IyBPAT8D0wMQyJvDSMrMZPT2ehWuSeO73Ha3oKj0Wi8hEEel7phlEpI6I3AOsBQYEL5opTTYfSOHmfy3mUGo6M0b0tqLLFEmR93ip6li30+p57qR7VDU+MLGMCayTGVmMnLqMxVsP87cbO3N7dHOvI5nAaQ88CXwmIjlAPLAXZ2DnOkBHoAMQCzyoqgu9CmrC18pdxxg6OZby5coxd3RfOjSq6XUkE6aK07n+ZWCzqr6pqm8CW0Xk74GLZkxgpKRlMmRSLEu2HuYfN3e1oquUUdVjqvonoAkwBlgH1AZaAVnAVKC7qp5nRZcpip82H2LQ+0uoXjmCD++xossUT3H6eF2uqo+dfqCqR0XkCuCJ/F4kIpWB74FK7vrnq+qzxchhzBmdyFTunBjLmj3J/PP2HlzTpZHXkUwJUdVTwHz3ZkxALEjcz/2zV9CqfjWmjYgmsmZlryOZMFecPl7lRKTG6QciUhOo4Mfr0oFLVLUr0A24UkT6FCOHMXk6nJrO2Ng01u09zrt39rSiywSMiNQWkfkisl5E1uXXv8yEr3lxu/jDzHg6NanJ3NF9rOgyAVGcPV5vAj+KyFz38a3A6wW9SFWVX86GrODe8roEkTFFduB4GndOXMq+EzlMGBrNRec08DqSCTJ3iJsrcA5BgtPva2GABlJ9E1igqgNFpCJQNQDLNCHk/e+38uLn67igbX3eG9yTqhWLdWljY/7Hr98kEemKMwZOErBIVTNUdZKIxAIXu7MNUtU1fi6vPE4H2DbAO6q6NI95RgGjACIjI4mJifFn0QCkpqYWan4vhVNWCI+8B0/m8HJcGsczlHs6Krp3DTF7vU7ln3DYvr5CNa+IjAD+BHyOU3AB9AaeFZFXVXViMZZdC7gQGAqgqhlARrECm5ChqryycAPjYrZwTZdGvH5LNypG2EhJJnAKLLzcAuhdnItgA2wSkUtUdY+qJgKJhV2pqmYD3USkNvAfEenkLst3nvHAeICoqCjt37+/38uPiYmhMPN7KZyyQujn3ZSUwuMTl5JBBHNG9yJ568qQzptbqG/f3EI472NAD1U94TtRRJ4GlgNFLrxwOu0fBCa7X0rjgQdyr8uEn+wc5amPEpkdu5NBvZvz1wGdKG/XXTQB5s8er8eAccBfgWbAG8BY4M7irlxVj4nIt8CVFKGAM8bXqt3HGDIplojy5Zg3ui/tzqpBzFavUxmPKFADyF0M1aD4XRsicI4A3KeqS0XkTZyTip4+PUNx9tgHWijslQyHDFk5ynur0onbn821Z1fg8tqH+OH774K2/mCwDKGRwZ/CqwXwqqoeAA6IyFBgdVFXKCINgEy36KoCXI5TyBlTZEu2Hmbk1GXUrlqBmSN706JeNa8jGW89CnwnIonAHndaU+Bc4JFiLns3sNuni8R8cp3NXZw99oEWCnslQz3DifQsxsyIJ27/SZ68ugN3X3h2UNcfLJYhNDL4U3iVB06dfqCqW0QEEWmkqvuKsM5GwFS3n1c5YJ6qflqE5RgDwDfrk7hnxnKa1a3KjBG9OauWnXlU1qnqpyLyBRANnL6Q3l4g1u3qUJxl7xeRXSLSTlU3AJfijIhvwtCxkxkMmxLHyl3HeHlgF26JauZ1JFPK+XuaxigR+RlIUNUjQDZQpSgrVNVVQPeivNaY3D5O2MMj81bSoVFNpg6Ppm61il5HMiHCLbAW554uIr3zOqGnkO4DZrpnNG4FhhVzecYDScfTuGtiLNsOnWDcHT25stNZXkcyZYA/hde3wMPA84CKyF6gMk4x9jWwLECnZxtTKDOX7uCpjxKJblmXCUOiqFHZn2HkjOEDoFiXL1DVBCAqIGmMJ7YfOsGdE5dy9EQGU4b1ol+b+l5HMmVEgYWXql4KICJnAz3dWw9gJE7HexWRraratiSDGuPr3ZgtjF2wnkvaN2TcHT2oXKG815FMCBGReWd6CqgbzCwm9Kzbd5zBE2PJzslh1t196NqstteRTBni94hwqroVZ5f6B6eniUhLnG99PQKezJg8qCovL9zAuzFb+H3Xxrx2S1cqlLcxdsxvXAYM5pfBmk8TnDG4TBm1bPsRhk2Jo3qlCOaM6kubhjUKfpExAVSsoXhVdTuwHbs2mgmC7BzlmY8TmbnUxtgxBYoBUlT1+9xPiMiq4McxoeDbDQe4Z0Y8jWtVYfrI3jSpXaSuysYUi10DwYSFtMxsHpqbwBeJ+xlzUWsev7IdIlZ0mbyp6o35PHd5MLOY0LBkbxYTvlxG+0Y1mDIsmvrVK3kdyZRRVniZkJeSlsnd05axZOsRnrqmAyMvCPwYO8aY0mv64u28tyqd6FZ2Io7xnhVeJqQdSElj6KQ4Nial8PqtXbmhe1OvI5kQJiKTzvCUAmnAZmCuqobJ1TtNcagq//xmM699tZHuDcszdXi0nYhjPGeFlwlZOw6fYPDEWA6mpDNhSBT92zX0OpIJfQ2AC4AcfrkMWSecTvXxwI3A8yJygTskhCmlcnKUv362lsk/befGHk24pv5RK7pMSLDTwUxIStyTzE3v/kxKWiaz7u5tRZfx10/AF0BTVb1QVS/EuVTQ58CXOJdA+wz4h3cRTUnLys7h0fkrmfzTdoad15JXB3a1E3FMyLA9Xibk/Lz5EKOmx1OzcgTTRvWmTcPqXkcy4eMB4BJVPXl6gqqeFJEXga9V9WURGQss8iyhKVFpmdncO2sFi9Yl8fDl53DfJW3sRBwTUqzwMiHl89X7eHBOAi3rV2Xq8Gga1bLTvU2hVMe5Huy6XNPPcp8DOI61faVSSlomI6cuI3b7Ef464FwG923pdSRjfsMaHxMypi/ZwTMfJ9KzeR0mDulFrap25pEptP8AE0XkMSDOndYLeBn4t/s4GtjoQTZTgg6lpjN0cizr96Xwxq3dGNCtideRjMmTFV7Gc6rK64s28dbXm7isQ0P+eXsPqlS0TrCmSMYArwEz+KV9ywImAY+6j9cBdwc/mikpe46dYvCEpexNPsX7d0VxcXvrE2pClxVexlOZ2Tn8+d+r+SB+Nzf3bMrfbuxMhF0CyBSR27drjIg8ArR2J29R1RM+8yR4kc2UjM0HUhg8MZbU9Cymj+hNr5Z2KU4T2qzwMp5JTc/iDzOX8/3Gg9x/aVseuqytdYI1AeEWWnZpoFJu1e5jDJkUS/ly5Zg7qi8dG9f0OpIxBbLCy3gi6XgawybHsSEphbE3debWXs29jmRKCRGJBP4IdMQZOHUtME5VkzwNZgLq582HuHvaMupUq8iMEb1pWb+a15GM8Ysd0zFBtzEphRvH/cz2wyeYMCTKii4TMCJyHs7o9IOAUzij1d8BbBKRvl5mM4GzcM1+hk6Oo0mdKnx4Tz8rukxYsT1eJqgWbznMqOnLqFyhPPNG96VTk1peRzKly6vAbGCMquYAiEg54F84g6b28zCbCYB5y3bxxIer6NqsNpOH9qJ21YpeRzKmUKzwMkHzccIe/vTBKprXq8qUYb1oWqeq15FM6dMNGHq66AJQ1RwReQ1Y4VkqExATftjKC5+t44K29fnXnT2pVsn+hZnwY7+1psSpKv/6bitjF6wnulVd3h8cZWN0mZKSDLQCNuSa3go4FvQ0JiBUlVe/3MA7327hms6NeO3WrlSKsCFnTHiywsuUqOwc5dlPEpmxZCfXdmnEP26xBtOUqDn8MoDqz+6084CxOIcgTZjJzlGe/jiRWUt3cnt0c164vpNdd9GENSu8TIk5kZ7FA3NWsGjdAUZfeDaPX9mectZgmpL1GCA4A6ZGuPczgHeBJzzMZYogIyuHh+Yl8Nmqffyhf2v+9Lt2NuSMCXtBL7xEpBkwDYjEOdV7vKq+GewcpmTtSz7FiCnLWL//OH+57lyG9GvpdSRTBqhqBvCAiPwfvx5A9WQ+LzMh6GRGFqOnx/PDpkP8+er2jLqwdcEvMiYMeLHHKwt4RFWXi0gNIF5EvlLVtR5kMSVg9e5kRk6LIzUti4lDetnlO0yJEpFP/JgHAFW9rsQDmWI7djKD4VPiSNh1jJdv6sItvZp5HcmYgAl64aWq+4B97v0UEVkHNMEZ5NCEuYVr9vPgnATqVK3A/Hv60aGRjSRtStxhrwOYwDlwPI3BE2PZdugE4+7owZWdGnkdyZiA8rSPl4i0BLoDS/N4bhQwCiAyMpKYmBi/l5uamlqo+b0UTlnhzHlVlQXbs5i3IYNWtcpxfw8hacNyknKfWxZkpWX7hqpQyKuqwzwNYAJmx+ET3DlxKUdSM5g8rBfntanvdSRjAs6zwktEqgMfAg+q6vHcz6vqeGA8QFRUlPbv39/vZcfExFCY+b0UTlkh77yZ2Tk883Eiczfs4prOzpmLlSuExpmLpWH7hrJwy2tC17p9x7lrUixZ2TnMursPXZvV9jqSMSXCk8JLRCrgFF0zVfXfXmQwgZF8MpM/zIrnp82H+ePFrXnk8nZ25qIxplCWbT/C8ClxVK0YwewxfWnTsIbXkYwpMV6c1SjARGCdqr4W7PWbwNlx+ATDp8Sx88hJXhnYhZujrAOsMaZwYjYcYMyMeBrVqsL0EdF2RQtT6nmxx+s8YDCwWkQS3Gl/VtXPPchiiih22xHGzIgnR5XpI3rT5+x6XkcyxoSZT1bu5eG5CZwTWYNpI6KpX72S15GMKXFenNX4I86ghiZMzY7dydMfJdK8blUmDu1Fq/rVvI5kjAkz05fs4JmPE+nVsi4ThkRRs7JdRsyUDTZyvfFbZnYO09em8/XO1Vx0TgPeur07tapYY2mM8Z+q8vY3m/nHVxu5tH1D3rmjR8icjGNMMFjhZfxy5EQGf5y5nMU7sxjlXv7HrpdmjCmMnBzlhc/WMemnbdzQvQkvD+xChfLlvI5lTFBZ4WUKtGF/CiOnxZF0PJ27O1fkz1d38DqSMSbMZGXn8PiHq/lw+W6G9mvJM9d2tDOgTZlkhZfJ15dr9vPQ3ASqVYpg7qg+JG9d6XUkY0yYychW7pm5nK/WJvHw5edw3yVt7GLXpsyywsvkybcfRtemtXhvcBRn1apMzFavkxljwklKWiavxaex/shJnh9wLnf1bel1JGM8ZYWX+Y1TGdk8On8ln63ax/XdGvP3m7pY51djTKEdTk1n6OQ4Nh3N4c3bujGgWxOvIxnjOSu8zK/sPHyS0TPiWb//OE9c1Z7RF55thwSMyYOIlAeWAXtU9Vqv84SaPcdOMXjiUvYcPcV93StZ0WWMywov8z/fbjjAA7NXICJMGtqLi9s19DqSMaHsAWAdUNPrIKFm84FUBk9cSmp6FjNG9ubE9lVeRzImZNh5vIacHOXNRZsYPiWOJnWq8t97z7eiy5h8iEhT4BpggtdZQs2q3ce45b3FZGYrc0b1oVfLul5HMiak2B6vMi75VCYPz03g6/UHuLF7E168oTNVKlp/LmMK8AbwGJDn1ZxFZBQwCiAyMpKYmJigBcstNTU1aOtfdzibN5enUb2i8KeoyhzcuIKYjcHNcCZeZ/B6/ZYhdDJY4VWGrd9/nDHT49l99BR/ue5c7urbwvpzGVMAEbkWOKCq8SLSP695VHU8MB4gKipK+/fPc7agiImJIRjrX7hmP68vWkHL+tWZNrw3Z9WqHPQM+fE6g9frtwyhk8EKrzLqk5V7eXz+KqpXjmDOqD5E2eEAY/x1HnCdiFwNVAZqisgMVb3T41ye+WDZLh7/cBVdmtZmyrBe1K5a0etIxoQsK7zKmMzsHMZ+sZ4JP24jqkUdxt3Rg4Y1Kxf8QmMMAKr6f8D/Abh7vB4ty0XXhB+28sJn6zi/TX3eG9yTapXs34ox+bG/kDJkX/Ip7pu1gmU7jjK0X0v+fHUHKkbY+RXGmMJTVf7x5Ube/nYzV3c+i9dv7UalCOsfakxBrPAqI2I2HODheStJz8y2gQyNCRBVjQFiPI4RdNk5yjMfJzJz6U5u69WMF2/oTHm77qIxfrHCq5TLys7hjUWbePvbzbQ/qwbv3NGD1g2qex3LGBOmMrJyeHheAp+u2seYi1rz+JXt7KQcYwrBCq9S7MDxNO6bvYKl245wS1RT/nJdJxsqwhhTZCczsrhnxnK+23iQ/7uqPaMvau11JGPCjhVepdRPmw/xwJwVpKZn8erNXRnYs6nXkYwxYSz5ZCbDp8axYudRxt7UmVt7Nfc6kjFhyQqvUiY7R/nnN5t48+tNtG5QnVl39+GcyDzHeDTGGL8cOJ7GXZNi2XrwBOPu6MGVnRp5HcmYsGWFVymSdDyNh+Ym8POWw9zYvQl/vb6TndptjCmWnYdPcufEpRxKTWfS0F6c37a+15GMCWv2X7mU+GptEo/NX0laZg5jb+rMLVHNrMOrMaZY1u8/zuCJsWRm5zDr7j50a1bb60jGhD1PCi8RmQScvuxGJy8ylBZpmdm89Pk6pi3eQcdGNXnr9u60aWhnLRpjiid+xxGGTY6jasUIZo3uS1vrsmBMQHi1x2sK8DYwzaP1lwobk1K4b9YKNiSlMOL8Vjx2ZTsbwNAYU2wxGw4wZkY8jWpVYdrwaJrVrep1JGNKDU8KL1X9XkRaerHu0kBVmbF0Jy98upYalSOYMqwX/ds19DqWMaYU+O/KvTw8L4G2DWswdXg0DWpU8jqSMaVKyPbxEpFRwCiAyMhIYmJi/H5tampqoeb3UmGzpmYokxLTWX4gm871yzOycwTsW0vMvrUlF9J3/WG0bcHylrRwy2vyN2PJDp7+OJFeLeoyYWgUNStX8DqSMaVOyBZeqjoeGA8QFRWl/fv39/u1MTExFGZ+LxUm6zfrk/jLh6s5djKHp67pwPDzWlEuyJfpCKdtC5a3pIVbXpM3VWVczBZeWbiBS9s35J07elC5gnVbMKYkhGzhZX6Rmp7Fi5+tZXbsLtqfVYOpw6Lp2Lim17GMMaVATo7y0ufrmPDjNm7o3oSXB3ahQvlyXscyptSywivExW47wiMfJLD76CnGXNSahy5vax3ojTEBkZWdwxP/Xs38+N0M7deSZ67tGPS96MaUNV4NJzEb6A/UF5HdwLOqOtGLLKEqPSub177cyPgfttKsTlXmje5Lr5Z1vY5ljCkl0jKzuW/2Cr5am8SDl7XlgUvb2th/xgSBV2c13u7FesPFmr3JPDx3JRuSUhjUuzlPXt3BRqA3xgRMSlomo6bFs3jrYZ77fUeGntfK60jGlBn23zyEZGTl8M63mxkXs5naVSsyeWgvLm5vw0QYYwLncGo6QyfHsXbfcd64tRvXd2/idSRjyhQrvELEyl3HeGz+KjYkpXB9t8Y8+/tzqVOtotexjDGlyN5jp7hz4lL2HD3F+3f15JL2kV5HMqbMscLLY+nZ7hlFP2ylYY3KTBoaZY2hMSbgthxMZfCEpaSkZTF9RG+iW1mfUWO8YIWXh5ZsPcwzP50i6eRWbo9uzv9d3d4GLDTGBNzq3ckMmRxLOYE5o/twbuNaXkcypsyywssDx9MyeXnBemYs2UmDKsKsu3vTr3V9r2MZY0qhdYezeefbJdSqUoEZI3vTqn41ryMZU6ZZ4RVEqsqnq/bx/KdrOZSazojzWxFdOcmKLmNMifhyzX7+EZ9Gq/rVmT6iN2fVqux1JGPKPCu8gmTboRM883EiP2w6ROcmtZg4JIouTWsTE3PA62jGmFJofvxuHv9wFS1qlGPe6L52so4xIcIKrxKWlpnNv77bwriYLVQqX47nB5zLHb1bUN5GhzbGlJAJP2zlhc/WcV6begxuecqKLmNCiBVeJeiHTQd5+qNEth8+ye+7NubpazrQsKbt6jfGlAxV5R9fbuTtbzdzVaezeOO2biz+8QevYxljfFjhVQK2HzrBC5+tY9G6JFrWq8r0EdFc0LaB17GMMaVYdo7y7CeJzFiyk1ujmvHSjZ1tz7oxIcgKrwBKScvk7W82M+mnbVQsX47HrmzH8PNaUbmCXdTaGFNyMrJyeHheAp+u2sfoi87miSvb23UXjQlRVngFQHaOMj9+F68s3MCh1AwG9mzKY79rZ4cVjTEl7lRGNmNmxPPdxoM8cVV7xlzU2utIxph8WOFVTD9uOsTfvljHmr3H6dmiDpOG9qJL09pexzLGlAHJJzMZPjWOFTuP8vcbO3NbdHOvIxljCmCFVxGt3p3M2AXr+XHzIZrUrsJbt3fn910a2e59Y0xQHDiexl2TYtl68ARvD+rB1Z0beR3JGOMHK7wKaduhE7z65QY+W7WPOlUr8PS1HbmzT3MqRVg/LmNMcOw8fJI7Jy7lUGo6k4b24vy2NgizMeHCCi8/HE/LZPGWw3y5JomPE/ZQMaIc91/ShrsvPJsadm1FY0wQrd9/nLsmxpKRncPMkb3p3ryO15GMMYVghVceUtOziNt+hCVbDrN462ES9ySTo1C9UgSDejfnvkva0qBGJa9jGmPKmPgdRxk2OZYqFcszb3Rfzoms4XUkY0whWeGFc1biyt3HiFl/gB82H2LV7mSyc5QK5YXuzepw7yVt6Xt2PXq2qEPFiHJexzXGlEHfbTzImOnxRNasxPQRvWlWt6rXkYwxRVBmC6/kU5l8v/Eg364/QMzGgxw5kUE5gS5NazP6wrPp27oeUS3qUqWi9d0yxnjr01V7eWhuAm0b1mDq8Gjb425MGCtThdeJ9CwWrUvivyv38t3Gg2RmK7WrVqD/OQ24uH1DLmzbwK5pZowJKTOX7uCpjxKJalGHCUN6UauK9Ss1JpyV+sJLVVm+8ygzluzki8R9pGXmcFbNygzp25IrO51F9+Z17LIaxpiQo6qMi9nCKws3cEn7hrwzqIftgTemFCjVhdfynUd58j+JrNt3nOqVIripR1Ou796Ens3rUM6KLWNMEYhIM2AaEAkoMF5V3wzkOlSVlz5fx/s/bGNAt8a8enNXKpS3/qXGlAaeFF4iciXwJlAemKCqfw/0OubE7uTpjxOJrFmZl27ozIBujalWqVTXmcaY4MgCHlHV5SJSA4gXka9UdW1AFp6dw//9ezUfxO9mSN8WPPv7c+2LojGlSNArEREpD7wDXA7sBuJE5JNANVppmdlMSkzn+92ruaBtff55e3dqV7V+W8aYwFDVfcA+936KiKwDmgDFbsPSMrO5f/YKvlybxAOXtuXBy9ra1TCMKWW82AUUDWxW1a0AIjIHGEAAGq247Ue4d9Zyko5n8Yf+rXnkinbWf8sYU2JEpCXQHViaa/ooYBRAZGQkMTExBS4rI1t5PT6NdUdyuKN9RbpX2Mt33+0tdsbU1FS/1l+SLIP367cMoZNBVDW4KxQZCFypqiPdx4OB3qp6b675fBuunnPmzClw2cnpyvur0rmscRbdmlQPfPgSkJqaSvXq4ZEVLG9Js7y/uPjii+NVNapEFh4AIlId+A54UVX/fab5oqKidNmyZQUuT1V58qNEerWsww3dmwYsZ0xMDP379w/Y8ixDeK7fMgQ3g4icsf0K2U5PqjoeGA9Ow+XvRhrwu9D4YP0VTlnB8pY0yxseRKQC8CEwM7+iq5DL5KUbOgdiUcaYEObFaTJ7gGY+j5u604wxJuSJ0+lqIrBOVV/zOo8xJrx4UXjFAW1FpJWIVARuAz7xIIcxxhTFecBg4BIRSXBvV3sdyhgTHoJ+qFFVs0TkXmAhznASk1R1TbBzGGNMUajqj4CdtWOMKRJP+nip6ufA516s2xhjjDHGKzYUsjHGGGNMkFjhZYwxxhgTJFZ4GWOMMcYEiRVexhhjjDFBEvSR64tCRA4COwrxkvrAoRKKE2jhlBUsb0mzvL9ooaoNSmjZQVOE9ivQQuF3yjJ4v37LENwMZ2y/wqLwKiwRWRbKlxrxFU5ZwfKWNMtrAi0UPiPL4P36LUPoZLBDjcYYY4wxQWKFlzHGGGNMkJTWwmu81wEKIZyyguUtaZbXBFoofEaWwfv1g2U4zdMMpbKPlzHGGGNMKCqte7yMMcYYY0KOFV7GGGOMMUEStoWXiFwpIhtEZLOIPJHH85VEZK77/FIRaelBTN88BeUdKiIHRSTBvY30IqebZZKIHBCRxDM8LyLylvteVolIj2BnzJWnoLz9RSTZZ9s+E+yMufI0E5FvRWStiKwRkQfymCdktrGfeUNqGxsQkXY+n0eCiBwXkQeDnOEh93cmUURmi0jlYK7fzfCAu/41wXr/ebVJIlJXRL4SkU3uzzoeZLjZ3Q45IlLiwymcIcMrIrLebdf+IyK1PcjwV3f9CSLypYg0LskMv6GqYXcDygNbgLOBisBKoGOuef4A/Mu9fxswN8TzDgXe9nrbulkuBHoAiWd4/mrgC0CAPsDSEM/bH/jU6+3qk6cR0MO9XwPYmMfvQ8hsYz/zhtQ2tttvPsPywH6cQR2Dtc4mwDagivt4HjA0yO+7E5AIVAUigEVAmyCs9zdtEvAy8IR7/wlgrAcZOgDtgBggyqPtcAUQ4d4f69F2qOlz//7TtUKwbuG6xysa2KyqW1U1A5gDDMg1zwBgqnt/PnCpiEgQM/ryJ2/IUNXvgSP5zDIAmKaOJUBtEWkUnHS/5UfekKKq+1R1uXs/BViH80/KV8hsYz/zmtB2KbBFVYM9gn4EUEVEInCKn71BXn8HnC8tJ1U1C/gOuLGkV3qGNsn3f9JU4PpgZ1DVdaq6oSTX60eGL93PAmAJ0NSDDMd9HlYDgnqWYbgWXk2AXT6Pd/PbfwT/m8f9kJOBekFJ91v+5AW4yd39OV9EmgUnWpH4+35CSV8RWSkiX4jIuV6HOc09BN4dWJrrqZDcxvnkhRDdxgZw9vrPDuYKVXUP8CqwE9gHJKvql8HMgLO36wIRqSciVXH2JHvVtkaq6j73/n4g0qMcoWQ4zp79oBORF0VkF3AHENSuEeFaeJVG/wVaqmoX4Ct++WZkim85ziGWrsA/gY+8jeMQkerAh8CDub6BhaQC8obkNjYgIhWB64APgrzeOjh7eVoBjYFqInJnMDOo6jqcw1lfAguABCA7mBnyos4xrjI9lpOIPAlkATO9WL+qPqmqzdz13xvMdYdr4bWHX39raepOy3Medzd3LeBwUNL9VoF5VfWwqqa7DycAPYOUrSj82f4hQ1WPq2qqe/9zoIKI1Pcyk4hUwCliZqrqv/OYJaS2cUF5Q3Ebm/+5CliuqklBXu9lwDZVPaiqmcC/gX5BzoCqTlTVnqp6IXAUp4+iF5JOdxdwfx7wKIfnRGQocC1wh1uEemkmcFMwVxiuhVcc0FZEWrnf5m4DPsk1zyfAEPf+QOAbDz/gAvPm6r9zHU4/mlD1CXCXe+ZdH5xDCPsKepFXROSs0/37RCQa5/feqyIcN8tEYJ2qvnaG2UJmG/uTN9S2sfmV2wnyYUbXTqCPiFR1fzcuxYN2TUQauj+b4/TvmhXsDC7f/0lDgI89yuEpEbkSeAy4TlVPepShrc/DAcD6YK4/IpgrCxRVzRKRe4GFOGfrTFLVNSLyPLBMVT/B+UcxXUQ243Ssuy3E894vItfh7Ho9gnOWoydEZDbOWWr1RWQ38CxQAUBV/wV8jtNXYjNwEhjmTVKHH3kHAveISBZwCrjN429Z5wGDgdUikuBO+zPQHEJyG/uTN9S2sQFEpBpwOTA62OtW1aUiMh/nMHQWsAJvLtXyoYjUAzKBP6rqsZJe4RnapL8D80RkBLADuMWDDEdwugI0AD4TkQRV/V2QM/wfUAn4yv2utkRVxwQ5w9Ui0g7IwfksSmz9eWayttEYY4wxJjjC9VCjMcYYY0zYscLLGGOMMSZIrPAyxhhjjAkSK7yMMcYYY4LECi9jTJHkdfHZfOa9UESWi0iWiAz0md7CnZ4gzsV7g3p2kTHGBJud1WiMKRIRuRBIxbmmZKcC5m0J1AQeBT5R1fnu9Io47VC6OzJ+ItBPVYN9TT9jjAkK2+NljCmSvC4+KyKtRWSBiMSLyA8i0t6dd7uqrsIZN8d3GRk+V2yohLVJxphSzho5EzJEpLeILBaRUyJyVESe9jqTKbTxwH2q2hNn79a4gl4gIs1EZBXORcHH2t4uUxJEZIqIfOp1DnCuYykiSSLSOkDL+0BEHgnEskzJs8LLhAQRuQz4DOeKA12Bl4HnRaSHp8GM39xDhf2AD9wR7t8DGuX7IkBVd7kXh28DDBGRyBINaoz3/gx8rqpbArS854EnRaRWgJZnSpAVXsZzbj+f94E/qeoEVd2oqn8D9gP9ReQaEXnb25TGD+WAY6razefWwd8Xu3u6EoELSiyhMR4TkarASJwvmQGhqquBrcCdgVqmKTlWeJlQcBFQG5iRa3omkA50ARKCG8kUlqoeB7aJyM3gXFxbRLrm9xoRaSoiVdz7dYDzgQ0lHtaUaSJSSUTecA/3pYnIEhE5P9c81URkmoikuvP9n4h8KiJTirn6qwEFfsojV3G6W3yCc0F0E+Ks8DKh4BJgtapmnp4gIg2BJjgX2O0CtHc7bK893WHbeMu9+OxioJ2I7HYv/nsHMEJEVgJrgAHuvL3cC9TeDLwnImvcxXQAlrrzfwe86n57N6YkvQzcCgwHugOrgQUi4nto/B84XwpvwGmjuhKYvbEXAPG5LyIfgO4WsUD06S8yJnTZcBLGcyKyAKijqr19pj2H80/8HJzDT+NV9U0RGQVEq+pIT8IaY8KSu6eqPk7BdRQYqarT3OfKAxuB2ar6lNtf8Qhwl6rOceepBuwGPlbVoe60/wD9ga9V1Xd8umtxCrdyOCeMTPB57iMgWVWH+EyriLOn93lVnewzfR/wiqq+JiLXAFep6r1neH9dgJVAmwD2HTMlwPZ4mVDQHWevyQgROUdE/gQ8DgwDKgJVgX+68ybgNJ7GGFMUrYEK+BzqU9VsnL23HXPNE+szzwmcL4G+3gTu8p0gIhHAazh7yboDfxKRej6zVAHSci2noO4WUHCXi1M+yzchzAov4ykRaQw0BAYB9+Hs8h8EDFDVH3EawnWqenr8px7AKi+yGmNKvUIdAlLVGCAl1+RoYI2q7lHVVOAL4Aqf5w8BdXK9pqDuFlBwl4u67s+DhXkPJvis8DJe64ZzJtzn7llwlVS1u6p+6T7fBWgtIhXchmgkv+z9MsaYwtoCZADnnZ7gHmrsC6z1mScT6OUzT1Ug3ys0uBoDe3we78EpoE5bwS971k7rjjOAsK8/4JypuMR93AXY5Y6R9wbOOHm+OgF7VDXJj4zGQxFeBzBlXnfy34PVBfgUiAPKAw+rqn2jM8YUiaqeEJF3gbEicgjYBjwEROIO+KuqqSIyyWeefcBTODsritsxeqG73Hqqetid1h2o5J6g8gPOSSmPA5erqopIJX7b5eLqXMu9wF22CXFWeBmv5Vt4qaqNxmyMCbTH3Z+TcfpWrQCuVNV9PvM8ClTDGaYhFXgdpzjL3T8rt738eg9XE37dV2y1iMQCtwHv+HS3uAZ4Caf4W8sv3S2ggC4XIlIZ5+zL3xX0xo337KxGY4wxpgDuXqcdOGcZ/sNnen/g3tNnNbqd69fhnO2YDMTjXPj9sM9rrsTpmN8Rp1iaqaq5+335rnsIzmj3nXD6h32Oc4bjQff5P+IUalecaRkmdNgeL2OMMSYXEemOM85cLFADZy9ZDWCuzzyLcMbcqnZ6nDpVXexeN/FbnEOTL/sWXQCqukBE3gGaUnB3Cyi4y0UmzslJJgzYHi9jjDEmF7fweh9oB2Th9Kt6VFXjA7ye+cA+VbXCqYywwssYY4wxJkhsOAljjDHGmCCxwssYY4wxJkis8DLGGGOMCRIrvIwxxhhjgsQKL2OMMcaYILHCyxhjjDEmSKzwMsYYY4wJkv8HhoWlMSUtp0AAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Set rho in [rho_ppoly_tab[0]/10,rho_ppoly_tab[2]*10]\n", "rho_for_plot = np.linspace(rho_ppoly_tab[0]/10,rho_ppoly_tab[NEOS-2]*10,1000)\n", "\n", "# Set P_{cold}(rho)\n", "P_for_plot = np.linspace(0,0,1000)\n", "for i in range(len(P_for_plot)):\n", " P_for_plot[i] = P_cold(rho_for_plot[i],rho_ppoly_tab,Gamma_ppoly_tab,K_ppoly_tab,NEOS)\n", "\n", "# Plot P_{cold}(rho) x rho\n", "import matplotlib.pyplot as plt\n", "f = plt.figure(figsize=(10,4))\n", "ax1 = f.add_subplot(121)\n", "ax1.set_title(r\"Plot #1: $P_{\\rm cold}\\times\\rho_{b}$\",fontsize='16')\n", "ax1.set_xlabel(r\"$\\rho_{b}$\",fontsize='14')\n", "ax1.set_ylabel(r\"$P_{\\rm cold}$\",fontsize='14')\n", "ax1.grid()\n", "ax1.plot(rho_for_plot,P_for_plot)\n", "\n", "ax2 = f.add_subplot(122)\n", "ax2.set_title(r\"Plot #2: $\\log_{10}\\left(P_{\\rm cold}\\right)\\times\\log_{10}\\left(\\rho_{b}\\right)$\",fontsize='16')\n", "ax2.set_xlabel(r\"$\\log_{10}\\left(\\rho_{b}\\right)$\",fontsize='14')\n", "ax2.set_ylabel(r\"$\\log_{10}\\left(P_{\\rm cold}\\right)$\",fontsize='14')\n", "ax2.grid()\n", "ax2.plot(np.log10(rho_for_plot),np.log10(P_for_plot))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 4: Taking a look at $\\epsilon_{\\rm cold}$ \\[Back to [top](#toc)\\]\n", "$$\\label{eps_cold}$$\n", "\n", "$\\epsilon_{\\rm cold}$ is determined via\n", "\n", "$$\n", "\\epsilon_{\\rm cold} = \\int d\\rho\\, \\frac{P_{\\rm cold}}{\\rho^{2}} \\ ,\n", "$$\n", "\n", "for some integration constant $C$. **Notation alert**: in the literature, the integration constants $C_{j}$ below are usually called $\\epsilon_{j}$. We will keep them as $C_{j}$ for a clearer exposition.\n", "\n", "In the case of a piecewise polytropic EOS, we then must have\n", "\n", "$$\n", "\\boxed{\n", "\\epsilon_{\\rm cold} = \n", "\\left\\{\n", "\\begin{matrix}\n", "\\frac{K_{0}\\rho^{\\Gamma_{0}-1}}{\\Gamma_{0}-1} + C_{0} & , & \\rho \\leq \\rho_{0}\\\\\n", "\\frac{K_{1}\\rho^{\\Gamma_{1}-1}}{\\Gamma_{1}-1} + C_{1} & , & \\rho_{0} \\leq \\rho \\leq \\rho_{1}\\\\\n", "\\vdots & & \\vdots\\\\\n", "\\frac{K_{j}\\rho^{\\Gamma_{j}-1}}{\\Gamma_{j}-1} + C_{j} & , & \\rho_{j-1} \\leq \\rho \\leq \\rho_{j}\\\\\n", "\\vdots & & \\vdots\\\\\n", "\\frac{K_{N-1}\\rho^{\\Gamma_{N-1}-1}}{\\Gamma_{N-1}-1} + C_{N-1} & , & \\rho_{N-2} \\leq \\rho \\leq \\rho_{N-1}\\\\\n", "\\frac{K_{N}\\rho^{\\Gamma_{N}-1}}{\\Gamma_{N}-1} + C_{N} & , & \\rho \\geq \\rho_{N-1}\n", "\\end{matrix}\n", "\\right.\n", "}\\ .\n", "$$\n", "\n", "We fix $C_{0}$ by demanding that $\\epsilon_{\\rm cold}\\left(\\rho=0\\right) = 0$. Then, continuity of $\\epsilon_{\\rm cold}$ imposes that\n", "\n", "$$\n", "\\frac{K_{0}\\rho_{0}^{\\Gamma_{0}-1}}{\\Gamma_{0}-1} = \\frac{K_{1}\\rho_{0}^{\\Gamma_{1}-1}}{\\Gamma_{1}-1} + C_{1}\\implies \\boxed{C_{1} = \\frac{K_{0}\\rho_{0}^{\\Gamma_{0}-1}}{\\Gamma_{0}-1} - \\frac{K_{1}\\rho_{0}^{\\Gamma_{1}-1}}{\\Gamma_{1}-1}}\\ ,\n", "$$\n", "\n", "for $C_{1}$ and\n", "\n", "$$\n", "\\boxed{C_{j} = C_{j-1} + \\frac{K_{j-1}\\rho_{j-1}^{\\Gamma_{j-1}-1}}{\\Gamma_{j-1}-1} - \\frac{K_{j}\\rho_{j-1}^{\\Gamma_{j}-1}}{\\Gamma_{j}-1}\\ ,\\ j\\geq1}\\ ,\n", "$$\n", "\n", "generically.\n", "\n", "\n", "\n", "## Step 4.a: *A function to evaluate $\\epsilon_{\\rm cold}$ based on `IllinoisGRMHD`* \\[Back to [top](#toc)\\]\n", "$$\\label{eps_cold__computing}$$\n", "\n", "We continue with our previous values of the piecewise polytrope EOS, for which we will have\n", "\n", "$$\n", "\\boxed{\n", "\\epsilon_{\\rm cold} = \n", "\\left\\{\n", "\\begin{matrix}\n", "\\frac{K_{0}\\rho^{\\Gamma_{0}-1}}{\\Gamma_{0}-1} & , & \\rho \\leq \\rho_{0}\\\\\n", "\\frac{K_{1}\\rho^{\\Gamma_{1}-1}}{\\Gamma_{1}-1} + C_{1} & , & \\rho_{0} \\leq \\rho \\leq \\rho_{1}\\\\\n", "\\frac{K_{2}\\rho^{\\Gamma_{2}-1}}{\\Gamma_{2}-1} + C_{2} & , & \\rho_{1} \\leq \\rho \\leq \\rho_{2}\\\\\n", "\\frac{K_{3}\\rho^{\\Gamma_{3}-1}}{\\Gamma_{3}-1} + C_{3} & , & \\rho \\geq \\rho_{2}\n", "\\end{matrix}\n", "\\right.\n", "}\\ ,\n", "$$\n", "\n", "remembering that the integration constants are set via\n", "\n", "$$\n", "\\boxed{\n", "C_{j} =\n", "\\left\\{\n", "\\begin{matrix}\n", "0 & , & {\\rm if}\\ j=0\\\\\n", "C_{j-1} + \\frac{K_{j-1}\\rho_{j-1}^{\\Gamma_{j-1}-1}}{\\Gamma_{j-1}-1} - \\frac{K_{j}\\rho_{j-1}^{\\Gamma_{j}-1}}{\\Gamma_{j}-1} & , & {\\rm otherwise}\n", "\\end{matrix}\n", "\\right.\n", "}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:06.155386Z", "iopub.status.busy": "2021-05-17T23:43:06.154479Z", "iopub.status.idle": "2021-05-17T23:43:06.156515Z", "shell.execute_reply": "2021-05-17T23:43:06.156958Z" } }, "outputs": [], "source": [ "# Setting eps_cold\n", "def eps_cold__integration_constants(rho_ppoly_tab,Gamma_ppoly_tab,K_ppoly_tab,NEOS):\n", " if NEOS == 1:\n", " eps_integ_consts[0] = 0\n", " return eps_integ_consts\n", "\n", " eps_integ_consts = [0 for i in range(NEOS)]\n", " for j in range(1,NEOS):\n", " aux_jm1 = K_ppoly_tab[j-1] * rho_ppoly_tab[j-1]**(Gamma_ppoly_tab[j-1]-1) / (Gamma_ppoly_tab[j-1]-1)\n", " aux_jp0 = K_ppoly_tab[j+0] * rho_ppoly_tab[j-1]**(Gamma_ppoly_tab[j+0]-1) / (Gamma_ppoly_tab[j+0]-1)\n", " eps_integ_consts[j] = eps_integ_consts[j-1] + aux_jm1 - aux_jp0\n", "\n", " return eps_integ_consts\n", "\n", "def eps_cold(rho,rho_ppoly_tab,Gamma_ppoly_tab,K_ppoly_tab,NEOS):\n", " if rho < rho_ppoly_tab[0] or NEOS==1:\n", " return K_ppoly_tab[0] * rho**(Gamma_ppoly_tab[0]-1) / (Gamma_ppoly_tab[0]-1)\n", "\n", " # Compute C_{j}\n", " eps_integ_consts = eps_cold__integration_constants(rho_ppoly_tab,Gamma_ppoly_tab,K_ppoly_tab,NEOS)\n", "\n", " # Compute eps_cold for rho >= rho_{0}\n", " if NEOS>2:\n", " for j in range(1,NEOS-1):\n", " if rho >= rho_ppoly_tab[j-1] and rho < rho_ppoly_tab[j]:\n", " return eps_integ_consts[j] + K_ppoly_tab[j] * rho**(Gamma_ppoly_tab[j]-1) / (Gamma_ppoly_tab[j]-1)\n", "\n", " if rho >= rho_ppoly_tab[NEOS-2]:\n", " return eps_integ_consts[NEOS-1] + K_ppoly_tab[NEOS-1] * rho**(Gamma_ppoly_tab[NEOS-1]-1) / (Gamma_ppoly_tab[NEOS-1]-1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 4.b: Plotting $\\epsilon_{\\rm cold}$ \\[Back to [top](#toc)\\]\n", "$$\\label{eps_cold__plotting}$$\n", "\n", "To create a reasonably interesting plot, which will also test all regions of our piecewise polytrope EOS, let us plot $\\epsilon_{\\rm cold}\\left(\\rho\\right)$ with $\\rho\\in\\left[\\frac{\\rho_{0}}{10},10\\rho_{2}\\right]$.\n", "\n", "**Remember**: We *don't* expect $\\epsilon_{\\rm cold}$ to be *smooth*, but we *do* expect it to be *continuous* (by construction)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:06.204126Z", "iopub.status.busy": "2021-05-17T23:43:06.186570Z", "iopub.status.idle": "2021-05-17T23:43:06.576031Z", "shell.execute_reply": "2021-05-17T23:43:06.576831Z" } }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAAEiCAYAAABX1AwVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABZjElEQVR4nO3dd5hU5fn/8ffNLm2pS+9FmlIEAUHUKMZGNBE1GlGj2KOJ35RfmsbEJCYmmt40iVFijQ0bltiiGzVK70Wk914Wdtlly9y/P85ZHNddmG1zZnY/r+uaiznnPOfM/czAmZtnnmLujoiIiIiIJE+jqAMQEREREWlolISLiIiIiCSZknARERERkSRTEi4iIiIikmRKwkVEREREkkxJuIiIiIhIkikJFxERERFJMiXhIiIiIiJJpiRcaszMrjIzj3vsN7MFZnazmWWGZX5iZlVeGcrMzjez/1eD2K4ws/Vx20vN7KvlyvQwsz+b2QdmdiCsQ5/qvqaI1D/14D53kZk9Y2brzKzAzJab2S/NrFUNXrda9U0mM/uTmb0UdRzxzOybZrbIzOo0B2vIn0+y3uOaSungJO1cDIwDvgjMBP4M3F7Da54PVPvLCRgFzAEws5bAoLLtOP2BLwF7gHdr8FoiUv+l633uO0Ap8ANgAvBX4CbgjVRPVKrLzPoBNwI/iTiU8v4OdAQmRx1IlOr480mL97he/sOTyMx39+nu/rq7Xw/kAN+IOKZDX07ASCAGLChX5h137+zu5wBPJzO4ZDGzrmZ2YiXHsszsc8mOSSRNpet97gvu/iV3f8zd/+vufwC+DowFxicr0CT7JrDA3WdHHUg8dy8AHib4j1GlGsB9+5vU0eeT6HscNSXhUpdmAa3NrFNlBcxsQtgNpMDMcs3seTMbFB57kOB/sd3jfgJem+iLh607I/j4y2k0sNTdC+PLuXusKpWqwusfZWYPmtkmMysys+1m9raZNa6L1zuCbwKvmdnJ5WJsDrwI/DNsQTssM8s0s++Y2eqwTjvM7Cd1ErFIekiX+9yOSmIH6J7o6yUQT6V1LVfuUjP70MwKw24D55lZjpnl1FIcTYEvA/+q4Fgq3JufAAZXlmSHvkkt3LfLnZvyn094vDa+axJ5jyOVGXUAUq/1Jfj5M6+ig2Y2AXgZeAu4BGgJ3AG8Z2YjgJ8R/Jx0PHBeeNrBI71o+AXWO27XK2YWf7ysj1xfd1+baGUs6Ce+Bvipu//kCGXbAv8DZhD85LsHaAd0d/fiRF+zFt0GHAX828wmuPv/4m7kw4HPunuFn1OZ8AtqGjCM4HNaD1wL/NjMlrn7k3VaA5HUlM73uVPDP5fFndeHBO9zFcR02Lq6+6aw3JnAYwT3k/9HUP8/AM2Aj6rymodxAtCWcl0MU+jePB/YT9A16P1KytT4vh0vHT6f8PVr67tmPkd+j6Pl7nroUaMHcBXgBP0QM4Fs4CsEX0zPh2V+Evx1+8R5s4EVQGbcvr5AMfC7cPtBYGMV4xlM0DL0O2BJ+HwEsA/4Vtx2kwrOvS6sS58KjvUGSoDbE4jhnPA6lwJZ8XWM8HPKBKaG78PpwJvATmB4guffAuQCveL2NQ6v96+o66eHHnX5qE/3ufD87sB24I1y+6tyn/tEfROpa7jvfWAxYHH7RoXvb04tfV7fJ+iW06Tc/pS5NxMkoK8foUy179vp+PmEx2rtuyaR9zjKh7qjSG36kOAf827gXoL/SV9TUUEza0HQd/FJdy8p2+/uawhaKU6t6LxEuPtSd58P9CS4YcwH8oFWwNPuPj98FFXxuuvcPdPd70ig+CyCm92/wtfeXpXXqkj4M3Xbwxzfa4eZ1SV8nycR3MTfJHj/z3D38n1HK7p2I4IWkXvc/dAsDB60HK0jaFEp66e418wyEqqUSPpJ+/tc2IXhBYJk++py163KfS7+mgnVNbw3jAae8TBLCsvNIWiBj7/mj8zsIzOLmdn55Y71M7P3wuPzzGx0uZC6AfsqqH+t35vDeCq9Px/m3rwjjLNSNblvl4shLT6fOviuOeJ7HCUl4VKbLiD4SfVooIW7X+nuuyspmw0YsKWCY1sJfh6sMjPLCPuSZQInAR+Ezz8DbAK2hsftsBequRbAvwl+7jwJOPnwxZMmk+AmFiNoWchK8LwxBD9Jvhy/M3wfu/PxzXkkwcC10lqJViT1pPV9Lq47w1HA2e6+sToxVCDRunYguPdUlPxuK7f9BkFXgncqKPs34CF3Hwh8D3isXH2bUXG3nlS6NxcAzRMoV937drx0+Xxq+7sm0fc4EkrCpTYtdvfZ7r7cyw0KqsAegp+2ulRwrAtBK1N1/IeglaoY6Ao8Ej5/gOAfcNmxardAHYmZtQb+C7zp7n9z9/fdfWm5MmPN7B0zWxg+Lo87dqaZzQ73/8/MRlbyOhPDgTMLzOzuBOJqBjxPcPM6geAG+qqZjUugWieEf5b/wj6R4Ob+fLg9BthsZi9aMA/xf8wsO4Hri6SLtL3PhX1tpxK0dJ7j7ouq+foVSbSuO8PYKhrI2jl+w4NZaFaXL2RmHQnuSQ+G5d4gSDBHxRXbRdDnOP68I96bw3IV3p/r4N7cjuD9qFQN79vxUv7zCdX2d80R3+MoKQmXSLh7PsFo/ovjf04ys94E/9hywl0Hqdr/Yr9C0Er1G2Bl+Px4gp+kfhi3XX4O3dp0CtALWF7RwfBG8QJBn8tjCQbYvBYe60gwovu68NgvgOfMrEm5a3QC/glc5O7DgdVAm8oCsmAk+nMEdT/D3WcRDMx5k+CGfkJl54bKbp794q7ZCPglsJDgi43w+j2AywhaCvcCn1g0RKShSKX7XPjv9THgs8D57j69WpWqRKJ1DVsuZwNfjG8ZNbNRBP2TE9EL2OKfHEi5Ntxf5kOgiZn1iNt32HtzGEeF9+c6ujf3PUIsNb1vH5Imnw/U/nfNYd/jqCkJlyj9CBgAvGRmXzCzSwn+l58L/DYssxRoZ2Y3mdnxZjbscBcMW6dmAwOBl8Pn+wl+YnsgbMGa7e7748+zYDW5i/j4BvC5cN+pcWV6m1mJmR1pYY6yqcAeNrMvmtl4M/uymU0J958IrHD3nDBmd/ey/6mfACwK+3fi7mU/yQ0s9xonAAvdfXG4/QDBALHK3EnQcnBG3LWLCW7o/wFetMNPdTWKoGXi3rBOF4bnHQdc6h9P8zgG+Lq77w/7E86iXOuJSAOTKve5ewgWGvotkG9mJ8Q9DiVCVbjPVbeuAD8GhhAkseeY2ZUEazRsJehyURvKukiMidt3pHszVH5/rtV7swX9xwdScVeOMjW9b5eX6p8P1OJ3TYLvcaSUhEtk3P1V4FyCn6SeIuhDtgw42d03h8XuJ2h9+AXB6nQvHum6YcvE6cCr4a7PAfPcfethTns6fNwYbt8bbv80/tJABkf4d+PuM4ArCPrV/RN4hWCVuk/9bJegRJcdPly5O4Hx7j7vEyd8fEM/zyuZ6sqCAT2Dwms8T7AS2SMEfe3Glf2ca2btCX4ajx8wdCJBq4pIg5RC97myhV1uAz4o97gu/tIkcJ+rSIJ1LeuecDlwDEFL7/eBbxMkebkJvNR6oKt9cl7vPuH+stdYS/BefiFuXyrdm88FigjqX5lq37crDCLFP586+K5J5D2OlqfAFC166NGQHgT92rYS3Fwh+NLrED7vSNB/7dhw+3MEP+M1IbiRt40rtwsYEm5fTyVTK9ZCvCeH1x51hHITCFpJjg63LwAWUckUaXrooYceZQ+CrgWFwI8qOJZD0IUmft9/gOvD52cSzHhi5cpcRZA0ZlUhjgrvz4e7N4fbTpDcJnRvJhgc+kjU73sqfT61/V2TDu+xhYGKSBKZ2ViCn/9ah7vucvd/hcfOBO4iGBGfB/yfu8+1YPGNbHffG5abCNxN0J/03wSt+CO8CgsQJRjrN4BfAa38MNM6mtmPCG7UAwkG9awBvuLhAhAiInBohpbf8fGc10cRzKDRmSB53RKW+wlBK31Hgu42hcAJ7r7RzAYADxEkyAeAG9x9ZrnXySRIzh5w999UIb4K78+V3ZvDcw7dn490b7ZgkaYZYV1XJhpXskT1+dTmd02qv8dllISLyGGZ2cMEN7JRRywsInIEYVeaJwn6T7cnmK/7XeAH/nFf6tp6rROAke5+b21etyYsWLky290fjzqWikT1+dTmd02qv8dllISLiIiIiCSZBmaKiIiIiCRZZtQBRKFDhw7ep0+fKp+Xn59PixYtaj+gCNW3OtW3+oDqlA6SWZ85c+bsdPeOSXmxFFHde3ZdS/W/x6kcn2KrvlSOT7FVrLL7doNMwvv06cPs2VWfNS0nJ4fx48fXfkARqm91qm/1AdUpHSSzPma2LikvlEKqe8+ua6n+9ziV41Ns1ZfK8Sm2ilV231Z3FBERERGRJFMSLiIiIiKSZErCRURERESSTEm4iIiIiEiSKQkXEREREUkyJeEiIiIiIkmmJFxEREREJMka5DzhIiLJVlwaY2tuIRv2HGDjngI27Sng7CFdGNytddShiYjIYWzfX8grC7fQPTuLMwd3rrXrKgkXEakFB0tK2by3kI17DrBpT0GQaO8tOLS9dV8hMf+4vBl0z26uJFxEJAXtPVDEq4u3Mm3BZqav3kXM4ZLRPZWEi4gkW2nM2ZJbwIbdBUFr9u4DbNhTwIbdB1i59QB7X331E+UzGhldWjejR3ZzTujXnh7ZWfRo25we2c3pnt2crm2a0yRTPQJFRFJF/sES3li6jRcXbOadFTsoLnX6tM/i5tP684Xh3RjQuVWtvp6ScBERwN3ZlV/EhrjkeuOeA4eS7s17Cygu/bgpu5FB1zZBUj20Qwajj+lLz+wsumcH+7q0bkZmhpJsEZFUtie/iDeXbeO1Jdt4d8UODpbE6NqmGVed2IfzhndnaPfWmFmdvLaScBFpMPIOlrBh9wHW7z4QJtkFYdIdPD9QVPqJ8u1bNKFHuyyGdW/DOcO60jM7i57tmtOrXdYnWrJzcnIYP35gFFUSEZEq2ry3gNeXbOW1JduYuXY3pTGnW5tmXDqmF+cM68ro3tk0alQ3iXc8JeEiUm/EYs7WfYWs23WA9bvzWbfrAOt2BUn2ht0H2HOg+BPlWzTJoGe7LHq3b8HJ/TvSs13zMNHOokd2c1o01S1SRCTdxWLO2txS/vyfFbyxbBsLN+YC0L9TS2489SjOHtKFYd3b1FmLd2WS+g1jZhOAPwIZwP3ufle5402Bh4FRwC7gEndfa2btganA8cCD7n5zBdeeBhzl7kPruBoiEqGikhib9hawdlc+68Mke92ufNaFLdxFJbFDZTMbGT2ym9OzXRZD41qyyxLt7KzGSb/piohI3cs9UMy7K3eQszx47Mw7CHzE8J5t+d6EQZw9pAv9OraMNMakJeFmlgHcA5wJbARmmdk0d18aV+xaYI+79zezScDdwCVAIfAjYGj4KH/tC4G8Oq6CiCTJgaKSQ63YhxLsXQdYuyufzXsLPjHLSPPGGfRun8VRHVrw2aM70atdFr3bZ9GnfQu6tlG/bBGRhqA05izdvI93VuwgZ/l25q7fS2nMadO8MZ8Z0IGuvosbzjuFjq2aRh3qIclsCR8DrHT31QBm9gQwEYhPwicCPwmfTwX+Ymbm7vnAe2bWv/xFzawl8P+AG4Cn6i58EalN+wuLWbvzAGt25bN2Z/7HLdu7D7Bj/8FPlG2b1Zje7bIY2SubC47rTq92WfTp0ILe7bLo2KqpWrNFRBoYd2fF9jzeX7mT91ftYvrqXewrLAFgaPfW3HRqP8YP6siInm3JzGhETk5OSiXgkNwkvDuwIW57IzC2sjLuXmJmuUB7YOdhrvsz4LfAgcO9uJndQJCo07lzZ3JycqoSOwB5eXnVOi+V1bc61bf6QHrXqTjmbD/gbM2PsS0/xtYDzrb8GFvyStn36uufKJvd1OiUZRzduhGndGlMp6xGdMoyOmU1okVjA0qA3OCRB/l5sHRtBJWqQDp/RiIi6cDdWbMznxlrdvP+ql18sGpX2MUEemQ3Z8LQLpzYrwMn9mtPp9bNIo42MWk96sjMRgD93P1bZtbncGXd/T7gPoDRo0f7+PHjq/x6wQwIVT8vldW3OtW3+kDq16k05mzeW8Dqnfms2ZHHmp35rA5btjft+WTXkfYtmtC3Qys6Ze1n3LB+9G3fgr4dW9C7XQuaN8mIrhI1lOqfkYhIujlYUsriTbnMXruH2ev2MHfdHnblFwHQsVVTTu7fnhP7dWBcv/b0bJcVcbTVk8wkfBPQM267R7ivojIbzSwTaEMwQLMy44DRZraWoC6dzCzH3cfXVtAi8vEc2qu2B0l2WaK9ZmfQhaSo9OPBkC2bZtK3QwtG9MzmguN6cFSHFvTt0II+HVrQpnljoCxp/VTvMhERaaB25xcxe+1u5qwLku5FG3MPfbf0aZ/FqYM6Mrp3O8b0zaZfx5b1ohtiMpPwWcAAM+tLkGxPAi4rV2YaMBn4ALgIeMvdnUq4+1+BvwKELeEvKQEXqb7SmLNxzwFWbs9j1Y688M98Vm7PI7fg4+n9mmQ0OjQY8vRjOgUt2h2CVu2OLdVHuyEys18DXwCKgFXA1e6+t4Jyh50lS0TqP3dn7a4DzFq7mzlr9zBr3W5W78gHoHGGMax7Gyaf2JtRvdsxqnd2yvXlri1JS8LDPt43A68R3HynuPsSM7sDmO3u04AHgEfMbCWwmyBRByBs7W4NNDGz84Gzys2sIiIJKigqZfXOj5PsVWHSvXpn/iem+OvQsgn9Orbk3GO70r9jS47q2IJ+HVvSrW1zMpKwkIGklTeAW8N7/d3ArcD34wskOEuWiNRDO/Yf5H8rd/LOih28t2In28MB+G2aN2Z072wuGtWD0b3bcWyPNjRrnL7dE6siqX3C3f0V4JVy+26Pe14IXFzJuX2OcO21VDB9oUhDtju/iJXb88q1bOexaW8BZb8xNTLo2S6L/h1bcsrAjvTv2JJ+nYJku21Wk2grIGnD3eNH2k4n+DWzvERmyRKReiAWc+Zt2MsbS7fxzkc7WLplHwDZWY05qX/Ql3tMn3b069gyKatTpqK0HpgpIoHcA8V8tH0/H23bz0db9/PRtjw+2rb/0CAWgGaNG3FUh5aM7JXNxaN60r9TkGz3ad+iwbQ6SNJcAzxZwf5EZskCamdGq7qW6rPipHJ8iq36Ujm+3H15/GXqf5i7rYS520vZe9DJMBiQ3YiLBjRmaIcMerVuRCPbBwX72LRsDZuWJSe2VHzflISLpJG8gyWs2LafFdvyWL4tTLq37Wfbvo/n1W7RJIMBnVtx+jGdGNi5VZBsd2xJ97bNG2xrg9QOM3sT6FLBodvc/YWwzG0E80k+VpPXqo0Zrepaqs+Kk8rxKbbqS7X43J2FG3N5du5GnpmdT15xIc0bZ3DqwM5MGNqF047udGhQfpRS7X0DJeEiKamwuJSV2/NYvnU//1lexMNrZ7F863427S04VKZZ40b079SSk/p3YGDnVgzq3IoBnYNkWwMjpS64+xmHO25mVwGfB06vZFB9IrNkiUga2LH/IE/P2cCzczexcnseTTIbMaJDBtecMYJTB3ZM62lnk0VJuEiE3J3NuYUs27yPD7fuY9mW/Szbso81u/IP9dnONOjfuYBRvbO5dExPBnZuxcDOrejZLkuDIyVlhLOefA841d0rWzwtkVmyRCRFuTvzN+zloffX8vKiLRSXOmP6tOOuC4fxuWFdmTfjf4wfWtGPZVIRJeEiSVJYXMpH24Ike9mW/Szdso8Pt+w7tMwuQK92WRzTtRVfGN6NQV2CZHvd4lmc/tlTIoxcJCF/AZoCb4S/xEx39xvNrBvBVITnVDZLVnQhi0giSmPOSws388B7a1i4MZeWTTO5fGxvrhjXm34dW0YdXtpSEi5Sy9ydrfsKDyXbwZ/7WLMz/9DqkVlNMji6S5BsH921NYO7tmJQl9a0bPrpf5Ib1dotacDdK1x9yd03A+fEbX9qliwRSU0lpTFemL+Ze95eyeqd+RzVsQV3TBzChSN7VPh9JVWjd1CkBtyddbsOsHhzLos37WPJ5lwWb8plz4GPF7bpkd2cY7q25txjuzG4ayuO7tKaXu2yNEhSRERSUizmvLBgE79/YwXrdx/g6C6tuPfykUwY0kXfXbVISbhIgkpjzpqd+SzZnMuijbks3pzLks372B92J2mcYQzq0oqzh3RhcLfWHNO1NYO6tKJ1s+hHhYuIiCRizrrd3PHiUhZszGVIt9bcd8Uozjims5LvOqAkXKQCJaUxVu3IZ9GmoGV7SZhwHygqBaBpZiOO6dqaiSO6MbRbG4Z2b8PAzq1oktko4shFRESqbtPeAn75yjJeWriFzq2b8tuLh3PBcd2VfNchJeHS4Lk7G3YXMH/jXhZs2Mv8DXtZsjmXwuJg+fasJhkM7tqaL43uydDubRjavTX9OrakcYYSbhERSW+lMefB99fy29eXE3Pn66cP4MZTjyKriVLEuqZ3WBqc3flFLNi4l/nr97IgTLzL+nA3zWzEsO5tuGxMb47tESTcfTu01FSAIiJS7yzZnMutzy5i4cZcThvUkZ+dP5Qe2VlRh9VgKAmXeq2wuJQlm3OZt34vCzbmsmDDXtbvDqYwNoOBnVpx5uDODO/ZlhE92zKwcyu1cIuISL1WUFTKH/7zEfe/u4bsrCb85bLjOHdYVy30lmRKwqVe2VMY4+WFW5izbg9z1u1myeZ9lITzAnZr04zhPdty2dhejOjZlqHd22iKJRERaVDeXbGD255bzPrdB5h0fE9u/dwxtMnSBAJRUAYiaas05ny4dR9z1+1h9ro9zF67J1zWfS5NMxsxvGdbrj/lKEb2ymZ4jzZ0at0s6pBFREQisSvvIHe+vIxn523iqA4tePz6ExjXr33UYTVoSsIlbeQdLGHe+iDZnrNuD/PW7yE/nK2kU6umjO6TzSldSrnk9OMZ3LW1ZioREZEGz915du4mfv7yUvIOlvD1z/bnq6f1p1njjKhDa/CUhEvK2l9YzOy1e5i+ehfT1+xm8aZcSmOOGRzdpTUXjOzO6N7tGNU7mx7ZzTEzcnJyGNGzbdShi4iIRG7drnxue24x763cychebbnri8cysHOrqMOSkJJwSRm5BcXMWrObGWt2MSNMumMeLIIzomdbbjq1H2P6tuO4Xm1ppQVwREREKlRcGuOB99bwhzc/IrNRI352/lAuH9NLc36nGCXhEpn9hcVMX707aOlevYulW/bhDk0yGzGiZ1tu/uwATujbjuN6ZdO8iX42ExEROZIFG/Zyy7OLWLZlH2cP6cxPzxtKlzYaE5WKkpqEm9kE4I9ABnC/u99V7nhT4GFgFLALuMTd15pZe2AqcDzwoLvfHJbPAp4G+gGlwIvufkuy6iNVU1waY/6Gvby3YifvrdzJ/A17KY05TTMbMbJXNt88fSBjj2rHiJ5t1VdNRESkCvIPlvDb1z/iwffX0LFVU/725VFMGNol6rDkMJKWhJtZBnAPcCawEZhlZtPcfWlcsWuBPe7e38wmAXcDlwCFwI+AoeEj3m/c/W0zawL8x8w+5+7/ruv6yJG5O6t25PHuip38b+VOpq/eTd7BEhoZDOsRdC85eUAHjuvVlqaZSrpFRESq460Pt/Gj55ewObeAL4/tzXcnDKK1um2mvGS2hI8BVrr7agAzewKYCMQn4ROBn4TPpwJ/MTNz93zgPTPrH39Bdz8AvB0+LzKzuUCPOq2FHFZuQTHvrthBzvIdvLdiJ1v3FQLQu30WE0d04zMDOjDuqA6ak1RERKSGtu8v5KcvLuXlhVsY2LklU28cx6je7aIOSxKUzCS8O7AhbnsjMLayMu5eYma5QHtg55EubmZtgS8QdHep6PgNwA0AnTt3Jicnp2rRA3l5edU6L5XVtE7uzuY8Z8GOEhbsKGXF3hgxhxaNYXD7DCb0bMKQ9hl0zDJgF+zcxbydy2st/vL0GaWH+lan+lYfEUlt7s5Tszdw58vLKCyO8e0zB/KVU/tpat40Uy8GZppZJvA48Keylvby3P0+4D6A0aNH+/jx46v8Ojk5OVTnvFRWnToVFJXy/qqdvL18O29/uCNcIAeO6dqam8Z35LRBnRjRsy2ZESz/rs8oPdS3OtW3+ohI6tqSF2PSfdOZsWY3Y/u24xcXDqNfx5ZRhyXVkMwkfBPQM267R7ivojIbw8S6DcEAzSO5D1jh7n+ohTilArvzi3hz2TZeX7KVd1fs5GBJjKwmGZzcvwM3f7Y/pw3qpNHXIiIidaS4NMZ976zm9+8XkNWkmLu/OIyLR/XUtINpLJlJ+CxggJn1JUi2JwGXlSszDZgMfABcBLzl7n64i5rZzwmS9etqPeIGbvPeAl5fspXXlmxj5trdlMac7m2bc+mYXpx+TCfG9G2nAZUiIiJ1bNHGXL73zEKWbdnH8V0yuOfaU+nUSg1f6S5pSXjYx/tm4DWCKQqnuPsSM7sDmO3u04AHgEfMbCWwmyBRB8DM1gKtgSZmdj5wFrAPuA34EJhrZgB/cff7k1Wv+mb1jjz+vXgrry3ZysKNuQAM6NSSm07tx4ShXRjSrTXh+ywiIiJ1qKColD+8+RH/eHc1HVs15b4rRtFkx4dKwOuJpPYJd/dXgFfK7bs97nkhcHEl5/ap5LLKCGto094CXlqwmRcXbmbxpn0AjOjZlu9POJqzh3TmKPU1ExERSaoPVu3i1mcXsnbXAS4d05NbPncMbZo3Jifnw6hDk1pSLwZmStXtzDvIK4u28Mj0Ala8+hYAw3u25YfnHsO5x3ala5vmEUcoIiLS8OwrLOaXr3zI4zPX07t9Fv+6fiwn9usQdVhSB5SENyCFxaW8tmQrU+ds5H8rdxJz6NHS+O7Zg/j8sV3p3b5F1CGKiIg0WG8s3cYPn1/Ejv0HueGUo/jWGQNp3kRjr+orJeH1nLszb8Nenp69kZcWbGb/wRJ6ZDfnq+P784Xh3djy4RzGj+9/5AuJiIhIndiZd5CfTFvCSwu3cHSXVvzjytEc26Nt1GFJHVMSXk9t31fIM3M3MXXOBlbtyKd54ww+N6wLF43qwQl92x+a0miLupaJiIhEwt15bt4m7nhpKQcOlvKds4JFdxpHsM6GJJ+S8HrE3Zm+ejePTl/Ha0u2UhJzRvfO5u4vHsU5w7rSqpmWihcREUkFG/cc4AfPLeadj3Ywqnc2d39xGP07tYo6LEkiJeH1wP7CYp6bt4lHPljHiu15tGnemKtP6sOlY3ppZhMREZEUEos5j0xfx92vBj9F//S8IVxxQm8tutMAKQlPY2t35vPAe2t4du5G8otKGda9Db+66FjOG96NZo01kENERCSVrNy+n+8/s4g56/Zw6sCO3HnBUHpkZ0UdlkRESXgamrd+D/e9s5pXl2ylcaNGfGF4N64c15vhPdtGHZqIiIiUU1wa4285q/jzWyvJaprB7740nAuO667F7xo4JeFpwt15e/l2/vbf1cxcs5vWzTK56dR+XHViHzq11spZIiIiqWjhxr18b+pCPty6n88f25WfnDeEDi2bRh2WpAAl4SnO3fnPsu384T8fsXjTPrq3bc6PPj+YS47vScum+vhERERSUUFRKb9/8yPuD5ec/8eVozlzcOeow5IUoiwuRZW1fP/hzRUs3JhLr3ZZ/PqiYzn/uO6aukhERCSFvb9qJ7c+u4h1uw5w6Zhe3HrO0bTWDGVSjpLwFDR3/R5+9tJS5q3fS892zfnVF4/lgpFKvkUkdZnZr4EvAEXAKuBqd99bQbm1wH6gFChx99FJDFOkThWXxrjjxaU8Mn0dfdpn8fj1JzCuX/uow5IUpSQ8hWzeW8Ddr37IC/M307FVU3554TAuGtVDybeIpIM3gFvdvcTM7gZuBb5fSdnT3H1n8kITqXsFRaV87V9zeevD7Vx3cl++fdYgLTkvh6UkPAXEYs6jM9Zx178/pDTm/N9n+3Pjqf1ooT7fIpIm3P31uM3pwEVRxSKSbHsPFHHtQ7OZt34Pv7hgGJeN7RV1SJIGlOVFbOOeA3z7qQXMWLObzwzowC8uGEbPdpozVETS2jXAk5Ucc+B1M3Pg7+5+X/LCEql9W3ILmDxlJmt3HuCey0byuWFdow5J0oSS8AjlLN/ON5+cT0mp86svHsvFo3tozlARSVlm9ibQpYJDt7n7C2GZ24AS4LFKLnOyu28ys07AG2b2obu/U8Fr3QDcANC5c2dycnJqowq1Ki8vLyXjKpPK8dWX2LbkxfjN7ELyi51vjWxG813LyclZnjLxJZtiqxol4RFwd/781kp+/+ZHDOrcir9+eRR9O7SIOiwRkcNy9zMOd9zMrgI+D5zu7l7JNTaFf243s+eAMcCnkvCwhfw+gNGjR/v48eNrFHtdyMnJIRXjKpPK8dWH2BZs2Mv/e3AWjTIbM/X6MQzt3qbug6N+vHdRSMXYlIQnWUlpjB8+v5gnZm3gguO684sLhmnghoikPTObAHwPONXdD1RSpgXQyN33h8/PAu5IYpgiteLdFTv4yiNzaN+yCY9cM5Y+akiTakjqtBtmNsHMlpvZSjO7pYLjTc3syfD4DDPrE+5vb2Zvm1memf2l3DmjzGxReM6fLIX7cxSVxPjqY3N5YtYGbj6tP7/70nAl4CJSX/wFaEXQxWS+mf0NwMy6mdkrYZnOwHtmtgCYCbzs7q9GE65I9by4YDPXPDiLXu2yeObGE5WAS7UlrSXczDKAe4AzgY3ALDOb5u5L44pdC+xx9/5mNgm4G7gEKAR+BAwNH/H+ClwPzABeASYA/67LulRHacz51pPzeX3pNn78hcFcfVLfqEMSEak17t6/kv2bgXPC56uB4cmMS6Q2PfT+Wn7y4hKO79OOf1w5mjbNtQCPVF8yW8LHACvdfbW7FwFPABPLlZkIPBQ+nwqcbmbm7vnu/h5BMn6ImXUFWrv79LD/4cPA+XVZiepwd257bhEvL9rCD889Rgm4iIhIGnF3fvfGR/x42hLOOKYzD18zRgm41Fgy+4R3BzbEbW8ExlZWJlzwIRdoD1S2qEP38Drx1+xeUcHaGGlf3ZG1r64p5onlRXzhqMb0L11PTs76Kl+jrqTiaOGaqG/1AdUpHdS3+ojIx0pjzo9eWMy/ZqznS6N78IsLhpGpRfSkFjSYgZm1MdK+OiNr31+5k6dem8GEIV3405dHptwUhKk4Wrgm6lt9QHVKB/WtPiISOFhSyjefmM+/F2/lq+P78d2zB6Xc97ikr2Qm4ZuAnnHbPcJ9FZXZaGaZQBtg1xGu2eMI14zMnvwivvHkfPp1bMlvvjRc/3BFJBJm1hToBjQHdrj7johDEkl5+wuL+cojc3h/1S5+eO4xXPeZo6IOSeqZZP6eMgsYYGZ9zawJMAmYVq7MNGBy+Pwi4K3K5poFcPctwD4zOyGcFeVK4IXaD716bp+2hD35Rfxh0ghaagl6EUkiM2tlZjeZ2TtALrASWAxsNbP1ZvYPMzs+2ihFUtOO/Qe59B/TmblmN7+/ZLgScKkTScsMwz7eNwOvARnAFHdfYmZ3ALPdfRrwAPCIma0EdhMk6gCY2VqgNdDEzM4HzgpnVvkq8CBBC8+/SZGZUV5bspUXF2zm22cOZEi35EzgLyICYGb/D7gNWE3QuHEnsBkoANoRzDL1GYLpBKcD/+fuKyIKVySl7DgQ4+K/vc/WfYX8Y/JoThvUKeqQpJ5KavOsu79CMI1g/L7b454XAhdXcm6fSvbP5tPTFkaqsLiUn7+8lIGdW3LT+H5RhyMiDc8JBIvmLK7k+ExgipndSDA17KmAknBp8JZt2cfPZxRCo0weu+4ERvXOjjokqcfUR6IOTPnfGjbsLuCx68ZqBLWIJJ27fynBcgeBe+s4HJG0MHPNbq59aBaNDZ68cRwDOreKOiSp55SE17L9hcX8/b+rOf3oTpzUv0PU4YiIiMgRvLF0Gzf/ay49spvz1cExJeCSFErCa9mj09eTW1DMN84YEHUoIiKHmNk5hzsedhcUaXCemrWBW55dyLAebfnnVcezcNb7UYckDYSS8FpUWFzK/e+u5tSBHTm2R9uowxERiVc23qYzMA74D2DAacAHlBuvI1LfuTt/++9q7n71Qz4zoAN/+/IoWmgmM0ki/W2rRS8v3MKu/CK+cqqmMhKR1OLuVwOY2ZvAMe6+NdzuAjwaZWwiyRaLOb94ZRn3v7eG84Z34zcXD6dJpsZwSXIpCa9Fj89cz1EdWjDuqPZRhyIiUpkewM647V18ctEzkXqtuDTG96cu5Nl5m7jqxD7c/vnBNGqkxfQk+ZSE15KPtu1n9ro93HbOMVoZU0RS2RPA/8zsOcCB84HHI41IJEkOFJXwtcfm8vbyHXznrIF87bT++s6WyCgJryXPzNlI4wzjwpHdow5FRKRS7v4TM3sZODHcdbO7z4kyJpFk2HugiGsenMX8DXv55YXDuHRMr6hDkgZOSXgtcHdeWbyFk/p3oH3LplGHIyJyWO4+C5gVdRwiybIlt4ArH5jJut0HuPfykUwY2jXqkESUhNeGJZv3sWF3ATef1j/qUEREKmRmswi6n3zqEODuPibJIYkkxcrteVz5wAz2FZbw0NVjGNdP47YkNSgJrwWvLdlKRiPjzMFdog5FRKQyF0UdgEiyzd+wl6v/OZOMRo144oYTGNq9TdQhiRyiJLwW/PejHYzs1ZZ2LZpEHYqISIXcfV3Z83BawuPDzZnuvi2aqETqzjsf7eDGR+fQoWVTHrl2DL3bt4g6JJFP0KSYNZR7oJhFm3I5sZ+WqBeR1GdmlwHvAecCnwfeNbNJ0UYlUrtemL+Jax+aRe/2LZh64zgl4JKS1BJeQx+s3oU7nNRfSbiIpIXvA8e7+x4AM8sGcgimLhRJew/+bw0/fWkpx/dpx/2TR9O6WeOoQxKpkJLwGvpg1U6aN85gRM+2UYciIpKIRkBe3HYe+lVU6gF35/dvfMSf3lrJWYM786dLj6NZ44yowxKplJLwGpq9bg+jemdruVsRSRePAu+b2TPh9oXAIxHGI1JjpTHnh88v5vGZ67lkdE/uvGAomRn6XpbUpiS8BgqLS1m+dT9fOfWoqEMREUmIu99tZv8BTgp33aTFeiSdFRaX8s0n5vPqkq187bR+fOesQVoFU9KCkvAa+HDrfkpizjBNeSQiacLMfgX8wt1nh9vZZnaXu98ScWgiVba/sJjrH57N9NW7uf3zg7nm5L5RhySSMP1WUwOLNu4FYFiPtpHGISJSBWe6+96yjXCA5lnRhSNSPTv2H2TSfdOZvXYPf7hkhBJwSTsJt4Sb2TmHO+7uryRwjQnAH4EM4H53v6vc8abAw8AoYBdwibuvDY/dClwLlAJfd/fXwv3fAq4jWAluEXC1uxcmWq+aWLgxl/YtmtCtTbNkvJyISG1oZGat3H0/gJm1BjR9hKSV9bsOcMWUGWzfd5D7J49m/KBOUYckUmVV6Y5ycfhnZ2Ac8B+C5Y5PAz4ADpuEm1kGcA9wJrARmGVm09x9aVyxa4E97t4/nLf2buASMxsMTAKGAN2AN81sINAF+Dow2N0LzOypsNyDVahXtS3alMvQ7m3U90xE0skfgffM7Mlw+xLg9xHGI1IlSzfvY/I/Z1JcGuOx68cysld21CGJVEvC3VHc/Wp3vxpoAhzj7he5+xeBwUDTBC4xBljp7qvdvYhgTtqJ5cpMBB4Kn08FTrcgw50IPOHuB919DbAyvB4E/5FobmaZQBawOdE61URpzFm9M59BXVol4+VERGqFu08BLgf2h4/Lwn0iKW/G6l1c8vcPyGxkTL1xnBJwSWvVGZjZA9gZt70r3Hck3YENcdsbgbGVlXH3EjPLBdqH+6eXO7e7u39gZr8B1gMFwOvu/npFL25mNwA3AHTu3JmcnJwEQv6kvLy8Q+dtPxCjqCRG8e6N5OSk74rP8XWqD+pbfUB1SgepXB8za+Tusfh97r4YWBxRSCLV8vqSrdz8+Dx6ZjfnkWvH0q1t86hDEqmR6iThTwD/M7PnCPphnw88XptBJSpc6W0i0BfYCzxtZl9290fLl3X3+4D7AEaPHu3jx4+v8uvl5ORQdt7bH26Hd2Zx7smjGN2nXXWrELn4OtUH9a0+oDqlgxSvz3tm9kt3fzF+p5k1d/eCqIISqYqnZm3glmcXcmyPtvzzquPJbtEk6pBEaqzKs6O4+0+AmwlanguBm939pwmcugnoGbfdI9xXYZmwe0kbgpb2ys49A1jj7jvcvRh4FjixilWqllU7ggXn+nVsmYyXExGpriHACgAzG2wfD2K5NG7BHpGU5O7cm7OS7z2zkJMHdORf149VAi71RrWmKHT3We7+x/CR6CIPs4ABZtbXzJoQDKCcVq7MNGBy+Pwi4C1393D/JDNramZ9gQHATIJuKCeYWVb4xXI6sKw6daqqVTvyadeiiW4GIpLqMoH88Pl0gl8OAWaQpEYLkeqIxZyfv7yMX726nIkjunH/laPJaqLlTaT+qMoUhbMIup986hDg7j6mgmOHhH28bwZeI5iicIq7LzGzO4DZ7j4NeAB4xMxWArsJEnXCck8BS4ES4GvuXgrMMLOpwNxw/zzCLid1bf3ufHq1y0rGS4mI1MRyYJyZ5QEtgbbh/ry45zVmZj8j6B4YA7YDV7n7pwbKm9lk4Ifh5s/d/aHyZUSKS2N8b+pCnpu3iatP6sOPzh1Mo0aaiUzql6r8l/Kimr5YOJf4K+X23R73vJCPp0Isf+6dwJ0V7P8x8OOaxlZVm/YUMEQrZYpI6vsTQePEemAOcD1wE3AqUJujyn/t7j8CMLOvA7cDN8YXMLN2BPfr0QSNOnPCqWr31GIckuYOFJXw1cfmkrN8B989exBfHd9PUwFLvZRwEu7u68qem1kX4Phwc6a7p+/0INUQizmb9xZy9pAuUYciInJY7v6gme0CBhH82vi4ma0HOhEk6LX1OvviNltQ8S+nZwNvuPtuADN7A5hARIP7JfXsyS/imodmsWDDXu66cBiTxvSKOiSROlPlzlVmdhlwB/AmQVeU35rZ7e7+RG0Hl6p25h2kqDRG92xNjyQiqS+cGeVFADM7F7iAYH2HWr1vm9mdwJVALsFCbuVVNFVt90quVeNpZetaKk9NCakdX0Wx7SqI8dvZhWwvcL42oildDqwmJ2d1SsSWSlI5PsVWNdUZ4fB94Piynw/DaQJzqOWbeSrbtDeY1au75igVkRRkZn3Dhc0+JRxPMzWurAE93H1DReXLXfdNgpWKy7vN3V9w99uA28zsVoJZtKrdVbA2ppWtayk+NWVKx1c+tpXb8/jBAzPYX5LBo9eN5oSj2qdMbKkmleNTbFVTndlRGhEM6CmTV83rpK1DSbhawkUkNX1gZg+Y2bjKCphZtpndRDDgvfzqxRVy9zPcfWgFjxfKFX0M+GIFl0hkqlppYOat38PFf3ufolLnia+cEGkCLpJM1WkJfxR4P25+2QuBR2ovpNS3aY9awkUkpR0N3Aa8bGYxggGZmwnWdsgGBgPHEEz1+k13f62mL2hmA9x9Rbg5EfiwgmKvAb8If0EFOAu4taavLenrvx/t4MZH5tCxVVMeuXYMvdu3iDokkaSpchLu7neb2X+Ak8JdN1VhrvB6YdPeAlo1y6RVs8ZRhyIi8inuvhf4rpndDpwLnAz0BpoDO4GHgNfC5etry11mNohgisJ1hDOjmNlo4EZ3v87dd4dTGc4Kz7mjbJCmNDwvzN/Et59awMDOrXjwmuPp1KpZ1CGJJFV1Bmb+CviFu88Ot7PN7C53v6XWo0tRm/YUqBVcRFJeuCz9VOL6gNfha1XU/YTwu+K6uO0pwJS6jkdS2xtri3nsw/mM7duOf0weTWs1akkDVJ2+3GeGrSwAhAM0z6q1iNLAltxCuikJFxERqRJ35zevLeexD4s4a3BnHrpmjBJwabCq0ye8kZm1cvf9AGbWGmhQ/4K27z/IsT20UI+IpJewL/ZZfDwt4GaCbilaLEfqXGnM+eHzi3h85gZO6ZHJvZePJDOjQc3rIPIJ1fnb/0fgPTP7gZn9AHgX+H3thpW6SmPO7vyDdGrVNOpQREQSZmbXAh8AYwnu/Y3C5++Hx0TqTGFxKV97bC6Pz9zAzaf15+ohTZSAS4N3xJZwM2vk7rGybXefYmYz+XghhsvcfUldBZhqduUdJObQUUm4iKSX7wEj3T0/fqeZ/QiYS7Capkit21dYzA0Pz2b66t38+AuDufqkvuTkbIk6LJHIJdId5T0z+2W44hoA7r7YzFaFg34alO37DwLQUaO4RSS9ONAKyC+3vxUVLzEvUmM79h9k8pSZfLRtP3+cNIKJIypcIFWkQUokCR8CrAAws8HAMnd34FIzO7eyEfH11Y5DSbhawkUkrXwH+K+ZLebjBXJ6ENzjvx1ZVFJvrd91gCumzGD7voPcP3k04wd1ijokkZSSSBKeycctJ9OBEcBqYAZwZ92ElbrKknD1CReRdOLuL5nZv4ExQLdw92ZgZriUvUitWbI5l8lTZlESi/Gv68dyXK/sI58k0sAkkoQvB8aZWR7QEmgb7s+Le95gbN9fCKglXETST5hsf1B+v5mNdfcZEYQk9dD01bu4/qHZtGqWyRM3jKN/p1ZRhySSkhIZmvwn4D7gvwRLH18f7j8V2FZHcaWsHfsP0qpZJs0aZ0QdiohIbXk66gCkfnhtyVaunDKTzm2aMfWmE5WAixzGEVvC3f1BM9sFDCIYPf+4ma0HOhEk6A3KrvwiOrRUK7iIpBcze6qyQ0C7ZMYi9dMTM9fzg+cWcWyPtvzzquPJbtEk6pBEUlpCi/WEM6O8CGBm5wIXAE2BJ+outNSUW1BMm+YNam0iEakfzgCuIOhKGM+AU5IfjtQX7s69Oav49WvLOXVgR/765ZFkNanOWoAiDUuV/5WEfQqnVufFzGwCwWI/GcD97n5XueNNgYeBUcAu4BJ3XxseuxW4FigFvu7ur4X72wL3A0MJptm6xt0/1eextuw5UERHtYSLSPrJAfa7+zvlD5jZwuSHI/VBLOb8/OVlTPnfGs4f0Y1fXzycxlqERyQhSfuvqpllAPcAZwIbgVlmNs3dl8YVuxbY4+79zWwScDdwSTg14iSCqbS6AW+a2cDwPwR/BF5194vMrAmQVZf12JNfzED1cRORNOPuFx7m2JnJjEXqh6KSGN+buoDn52/mmpP68sNzj6FRI4s6LJG0kcz/ro4BVrr7ancvIujKMrFcmYnAQ+HzqcDpZmbh/ifc/aC7rwFWAmPMrA3Bz6gPALh7kbvvrctK5BYU0yZL3VFERKThOlBUwvUPz+b5+Zv53oRB/OjzSsBFqiqZnba6AxvitjcCYysr4+4lZpYLtA/3Ty93bnegANgB/NPMhhPM3vKN8ssyA5jZDcANAJ07dyYnJ6fKFdi7L4+8g8bebZvIydlR5fNTUV5eXrXei1RV3+oDqlM6SOX6mNmUSg45UEjQqPGku29OXlSSzvbkF3H1g7NYuHEvd104jEljekUdkkhaSveRE5nASOD/3H2Gmf0RuAX4UfmC7n4fwVSLjB492sePH1/lF3vhtbeBA4wcMpDx4/rUIOzUkZOTQ3Xei1RV3+oDqlM6SPH6dAQ+A8SAxeG+oQQDMucAFwJ3mNln3H1+JBFK2ti8t4Arp8xk/e4D/PXLozh7SJeoQxJJW8nsjrIJ6Bm33YOPl07+VBkzywTaEAzQrOzcjcDGuEUmphIk5XUir9gBaJOlaZdEJG38D/g30MPdT3H3Uwjuoa8ArwO9gZeB30YXoqSDldv388W/vs+23EIevmaMEnCRGkpmEj4LGGBmfcMBlJOAaeXKTAMmh88vAt5ydw/3TzKzpmbWFxhAsNTyVmCDmQ0KzzkdWEodySsKkvBs9QkXkfTxDeAOdz9QtiN8fifwrXCMzt3AiGjCk3Qwb/0eLvrbBxSXOk985QROOKp91CGJpL2kdUcJ+3jfDLxGMEXhFHdfYmZ3ALPdfRrBAMtHzGwlsJsgUScs9xRBgl0CfC2cGQXg/4DHwsR+NXB1XdUhv7gsCVdLuIikjZZAV2BZuf1dwmMA+0j/7olSR/770Q5ufGQOnVo35ZFrxtKrfZ1OQibSYCT1puvurxD8BBq/7/a454XAxZWceydBy035/fOB0bUaaCXKuqO0VUu4iKSP54AHzOx7BL9IAhwP/Ap4NtweA3wUQWyS4l6Yv4lvP7WAgZ1b8dA1Y+jYSutkiNQWtXxUQX5x8GdrrZgpIunjRuB3wKN8fM8vAaYA3wm3lwHXJz80SWX//N8afvriUk44qh33XTma1s303SdSm5SEV0FhSdAS3lLL8YpImgj7f99oZt8G+oW7V8VP5apZUSSeu/Ob15dzz9urmDCkC3+YNIJmjTOiDkuk3lE2WQWFpU5WkwwtSCAiaSdMurU8vRxWSWmMHz6/mCdmbeDSMb34+flDydB3nkidUBJeBQUl0KKp3jIRSS9m1hn4GjCYYJGepcC97r4t0sAkpRQWl/L1x+fx+tJt/N9n+/P/zhxIsGi1iNSFZE5RmPYKS5yWSsJFJI2Y2UkEq2JeRrDKcCFwObDCzMZFGZukjn2FxUyeMpPXl27jJ18YzLfPGqQEXKSOKaOsgsJSaNFU/eJEJK38BngcuNHdYwBm1gj4G8ECPSdGGJukgO37C5k8ZRYrtu3nj5NGMHFE96hDEmkQlIRXQWGJ066V3jIRSSsjgKvKEnAAd4+Z2e+AeZFFJSlh3a58rnhgJjvzDvLAVcdz6sCOUYck0mAoo6yCwhLUHUVE0k0u0BdYXm5/X2Bv0qORlLFkcy6Tp8yiNBbjsevGclyv7KhDEmlQlFFWQWGpa2CmiKSbJ/h4sZ73w30nESxV/3hkUUmkPli1ixsenk2rZpk8fMOJ9O/U8sgniUitUkZZBYUlSsJFJO18DzCCxXkyw+dFwF+BWyKMSyLy6uKtfP2JefRql8Uj146ha5vmUYck0iApo6yCAnVHEZE04+5FwDfM7FY+uVjPgQjDkog8MXM9P3huEcN7tmXK5OPJbtEk6pBEGixllAmKxZziGFo1TERSnplNS6AMAO5+Xp0HJJFzd+7NWcWvX1vO+EEduffykWRp9WeRSOlfYIKKSoOJBZo11tTqIpLydkUdgKSOWMy546WlPPj+Wi44rju/uuhYGmfou0wkakrCE3SwJEjCm+jGJSIpzt2vTvZrmtnPgIlADNhOMC3i5grKlQKLws31aomvW0UlMb7z9AKmLdjMtSf35bZzjqGRlqEXSQlKwhN0sKQUgKbqjiIiUpFfu/uPAMzs68DtwI0VlCtw9xHJDKyhOljiXPfwbN75aAffn3A0N556lFbBFEkhSsITVBS2hDfNVEu4iEh57r4vbrMF4FHFIrAnv4i7ZxWydt8B7v7iMC45vlfUIYlIOUrCE3RQSbiIyGGZ2Z3AlQQLBJ1WSbFmZjYbKAHucvfnK7nWDcANAJ07dyYnJ6fW462pvLy8lIxrV0GM38wuZMeBGDePaEbn/NXk5KyOOqxPSNX3DlI7Nkjt+BRb1SgJT9DBYiXhItKwmdmbQJcKDt3m7i+4+23AbeF0iDcDP66gbG9332RmRwFvmdkid19VvpC73wfcBzB69GgfP358rdWjtuTk5JBqca3Ytp9bp8wkrzSD7x7flK9ceHrUIVUoFd+7MqkcG6R2fIqtapKaUZrZBDNbbmYrzexTi0SYWVMzezI8PsPM+sQduzXcv9zMzi53XoaZzTOzl+oq9rLZUZpmqk+4iDRM7n6Guw+t4PFCuaKPAV+s5Bqbwj9XAznAcXUadAMyd/0eLv77B5TEnCdvGMegdvq+EkllSUvCzSwDuAf4HDAYuNTMBpcrdi2wx937A78nWFaZsNwkYAgwAbg3vF6ZbwDL6jL+g8XhwEy1hIuIfIqZDYjbnAh8WEGZbDNrGj7vAJwELE1OhPXb28u3c/k/ZtC2eWOeufFEBndrHXVIInIEycwoxwAr3X11uILbEwQ36ngTgYfC51OB0y0Yyj0ReMLdD7r7GmBleD3MrAdwLnB/XQZ/aIpCJeEiIhW5y8wWm9lC4CyCxhHMbLSZld2fjwFmm9kC4G2CPuFKwmvo+XmbuP6h2RzVsQVP33givdpnRR2SiCQgmX3CuwMb4rY3AmMrK+PuJWaWC7QP908vd2738PkfgO8BrQ734jUd5DNvewkAixfMI29t/fmJLxUHKtREfasPqE7poL7VpzrcvbLuJ7OB68Ln7wPDkhlXfffAe2v42UtLGXdUe+67chStmjWOOiQRSVBaD8w0s88D2919jpmNP1zZmg7yyV+4BebO5cQTxjCw82Hz/bSSigMVaqK+1QdUp3RQ3+ojqc/d+fVry7k3ZxUThnThD5NG0EzrWIiklWT2rdgE9Izb7hHuq7CMmWUCbQiWX67s3JOA88xsLUH3ls+a2aN1EXzZYj1aMVNERKJUUhrjlmcWcW/OKi4b24t7Lh+pBFwkDSUzo5wFDDCzvmbWhGCg5bRyZaYBk8PnFwFvubuH+yeFs6f0BQYAM939Vnfv4e59wuu95e5frovgD80T3lhJuIiIRKOwuJSvPjaXJ2dv4Ouf7c+d5w8lQ8vQi6SlpHVHCft43wy8BmQAU9x9iZndAcx292nAA8AjZrYS2E2QWBOWe4pgFH0J8DV3L01W7BC/YqZaG0REJPn2FRZz3UOzmbV2Nz89bwiTT+wTdUgiUgNJ7RPu7q8Ar5Tbd3vc80Lg4krOvRO48zDXziGYc7ZOHOqOotlRREQkybbvK2TyP2excvt+/jjpOM4b3i3qkESkhtJ6YGYyacVMERGJwrpd+VzxwEx25h3kgcnHc8rAjlGHJCK1QEl4gopKYxiQqb53IiKSJIs35XLVP2dSGnP+df0JjOjZNuqQRKSWKAlPUFFJjMxGEKwdJCIiUrc+WLWL6x+eTZvmjXnomjH079Qy6pBEpBYpCU9QzB01gouISDK8ungLX398Pr3bZ/HwtWPo2qZ51CGJSC1TEp6gmINycBERqWuPz1zPbc8tYkTPtky56njaZjWJOiQRqQNKwhPkDuqJIiIidcXdueftlfzm9Y84bVBH7r18FM2baFpckfpKSXiCYu5qCRcRkToRizl3vLSUB99fy4XHdefui46lsVZoFqnXlIQnyN3VEi4iIrWuqCTGd55ewLQFm7nu5L784JxjaKRBSCL1npLwBDnqEy4iIrUr/2AJNz46h3dX7OSWzx3NV045SrNwiTQQSsITpO4oIiJSm3bnF3H1g7NYtHEvv7roWL40umfUIYlIEikJT5AGZoqISG3ZtLeAKx6YwaY9Bfz9itGcObhz1CGJSJIpCU9QzEEdUkREpKZWbNvPFQ/MJL+ohEeuHcuYvu2iDklEIqAkPGFarEdERGpmzro9XPPgLJpkNuKpr4zjmK6tow5JRCKiJDxBsVjUEYiISDp7e/l2bnp0Dl1aN+ORa8fSs11W1CGJSISUhCdIy9aLiEh1PTdvI999eiFHd23Fg1ePoUPLplGHJCIRUxKeIE1RKCIi1XH/u6v5+cvLOLFfe/5+xShaNWscdUgikgKUhCco5h51CCIikkbcnV+9tpy/5qzinGFd+P0lI2iaqWXoRSSgJDxRjrqjiIhIQkpKY/zguUU8NXsjl4/txR0Th5KhLxERidMomS9mZhPMbLmZrTSzWyo43tTMngyPzzCzPnHHbg33Lzezs8N9Pc3sbTNbamZLzOwbdRW7WsJFRCQRhcWl3PjoXJ6avZFvnD6An5+vBFxEPi1pLeFmlgHcA5wJbARmmdk0d18aV+xaYI+79zezScDdwCVmNhiYBAwBugFvmtlAoAT4trvPNbNWwBwze6PcNWuFo5ZwERE5vNyCYq5/aDaz1u3mp+cNYfKJfaIOSURSVDJbwscAK919tbsXAU8AE8uVmQg8FD6fCpxuZhbuf8LdD7r7GmAlMMbdt7j7XAB33w8sA7rXRfAxNYSLiMhhbN9XyCV//4B5G/bwp0nHKQEXkcNKZp/w7sCGuO2NwNjKyrh7iZnlAu3D/dPLnfuJZDvsunIcMKOiFzezG4AbADp37kxOTk6Vgt+2rRA8VuXzUl1eXl69qlN9qw+oTumgvtVHqm7tznyumDKDXXlFTLnqeD4zoGPUIYlIiqsXAzPNrCXwDPBNd99XURl3vw+4D2D06NE+fvz4Kr3G05vmsmH/Vqp6XqrLycmpV3Wqb/UB1Skd1Lf6SNUs3pTLVf+cSczh8etPYHjPtlGHJCJpIJndUTYBPeO2e4T7KixjZplAG2DX4c41s8YECfhj7v5snUROMDDT1CdcRETivL9qJ5Pum07TzAyevnGcEnARSVgyk/BZwAAz62tmTQgGWk4rV2YaMDl8fhHwlrt7uH9SOHtKX2AAMDPsL/4AsMzdf1eXwbsneSoZERFJabO2lnDVlFl0a9uMZ246kX4dW0YdkoikkaR1Rwn7eN8MvAZkAFPcfYmZ3QHMdvdpBAn1I2a2EthNkKgTlnsKWEowI8rX3L3UzE4GrgAWmdn88KV+4O6v1Hb8mqJQRETKPDZjHffOP8jI3tk8MHk0bbOaRB2SiKSZpPYJD5PjV8rtuz3ueSFwcSXn3gncWW7feyRpNXkHTP1RREQaNHfnz2+t5HdvfMTwjhk8eu1YmjfRKpgiUnX1YmBmMrh7crJ9ERFJSbGY89MXl/DQB+u4cGR3zumwRwm4iFSbujknyB0NzBQROQIz+7aZuZl1qOT4ZDNbET4mV1QmFRWVxPjGk/N56IN1XP+ZvvzmouFkagU3EakBtYQnKKaWcBGRwzKznsBZwPpKjrcDfgyMJujlNydcOXlP8qKsuvyDJdz46BzeXbGTWz93NF85tV/UIYlIPaCW8AQFfcKjjkJEJKX9HvgewS2zImcDb7j77jDxfgOYkKzgqmN3fhGX/WM676/axa8uOlYJuIjUGrWEJyjmSRoBKiKShsxsIrDJ3RccZhB7RSsnd6+oYE1XOa4NOwti/GZ2IbsKnJtHNKVT3ipyclYdOp7qK6WmcnyKrfpSOT7FVjVKwhOkgZki0tCZ2ZtAlwoO3Qb8gKArSq2o6SrHNfXRtv3c8sBM8ksz+NcNx3N8n3afKpPqK6WmcnyKrfpSOT7FVjVKwhOkgZki0tC5+xkV7TezYUBfoKwVvAcw18zGuPvWuKKbgPFx2z2AnDoJtgbmrNvNNQ/OpmlmI56+cRxHd2kddUgiUg+pT3iCNDBTRKRi7r7I3Tu5ex9370PQzWRkuQQcgsXazjKzbDPLJmg5fy3J4R7W2x9u5/L7Z9CuRROeuelEJeAiUmeUhCdILeEiIlVnZqPN7H4Ad98N/AyYFT7uCPelhGfnbuS6h2fTv1NLnr5xHD3bZUUdkojUY+qOkiC1hIuIJCZsDS97Phu4Lm57CjAlgrAO6/53V/Pzl5dxYr/2/P2KUbRq1jjqkESknlMSniBNUSgiUv+4O3e/upy//XcV5wzrwu8vGUHTTK2CKSJ1T0l4gjQ7iohI/VJSGuMHzy3iqdkb+fIJvfjpeUPJ0CqYIpIkSsITpD7hIiL1R2FxKTf/ax5vLtvGN04fwDfPGMBh5jcXEal1SsITpD7hIiL1Q25BMdc/NJtZ63bzs4lDuGJcn6hDEpEGSEl4goIVM5WGi4iks+37CrlyykxW7cjjz5cex+eP7RZ1SCLSQCkJT5AGZoqIpLc1O/O54oEZ7M4v4p9XjeHkAR2iDklEGjAl4QnSwEwRkfS1eFMuk6fMxIEnbjiBY3u0jTokEWnglIQnyB00aF5EJP28v3InNzwyhzbNG/PItWM4qmPLqEMSEVESnqiYu5YXFRFJM68s2sI3n5hP3w4teOiaMXRp0yzqkEREgCQvW29mE8xsuZmtNLNbKjje1MyeDI/PMLM+ccduDfcvN7OzE71mbVFLuIhIenl0+jq+9q+5HNujDU99ZZwScBFJKUlLws0sA7gH+BwwGLjUzAaXK3YtsMfd+wO/B+4Ozx0MTAKGABOAe80sI8Fr1oqYe11cVkRE6sA9b6/kh88v5rODOvHItWNpk6Vl6EUktSSzJXwMsNLdV7t7EfAEMLFcmYnAQ+HzqcDpFqyeMBF4wt0PuvsaYGV4vUSuWStOGdiRAdnqkCIikg76dmjBl0b34G9XjKJ5Ey1DLyKpJ5l9wrsDG+K2NwJjKyvj7iVmlgu0D/dPL3du9/D5ka4JgJndANwA0LlzZ3JycqoU/IlZkNexqMrnpbq8vLx6Vaf6Vh9QndJBfatPfXDOsK6cM6xr1GGIiFSqwQzMdPf7gPsARo8e7ePHj6/yNXJycqjOeamsvtWpvtUHVKd0UN/qIyIidS+Z/Ss2AT3jtnuE+yosY2aZQBtg12HOTeSaIiIiIiIpJZlJ+CxggJn1NbMmBAMtp5UrMw2YHD6/CHjL3T3cPymcPaUvMACYmeA1RURERERSStK6o4R9vG8GXgMygCnuvsTM7gBmu/s04AHgETNbCewmSKoJyz0FLAVKgK+5eylARddMVp1ERERERKojqX3C3f0V4JVy+26Pe14IXFzJuXcCdyZyTRERERGRVKY590REREREkkxJuIiIiIhIkikJFxERERFJMvMGuBy7me0A1lXj1A7AzloOJ2r1rU71rT6gOqWDZNant7t3TNJrpYQa3LPrWqr/PU7l+BRb9aVyfIqtYhXetxtkEl5dZjbb3UdHHUdtqm91qm/1AdUpHdS3+khiUv1zT+X4FFv1pXJ8iq1q1B1FRERERCTJlISLiIiIiCSZkvCquS/qAOpAfatTfasPqE7poL7VRxKT6p97Ksen2KovleNTbFWgPuEiIiIiIkmmlnARERERkSRTEi4iIiIikmRKwitgZhPMbLmZrTSzWyo43tTMngyPzzCzPhGEmbAE6nOVme0ws/nh47oo4kyUmU0xs+1mtriS42Zmfwrru9DMRiY7xqpKoE7jzSw37jO6PdkxVoWZ9TSzt81sqZktMbNvVFAmrT6nBOuUVp+TVI+ZDYr7jOeb2T4z+2bUcZUxs2+Ff0cXm9njZtYs6pjKmNk3wriWpMJ7VtG918zamdkbZrYi/DM7hWK7OHzvYmYW6XR7lcT3azP7MLynP2dmbVMotp+Fcc03s9fNrFsUsX2Cu+sR9wAygFXAUUATYAEwuFyZrwJ/C59PAp6MOu4a1ucq4C9Rx1qFOp0CjAQWV3L8HODfgAEnADOijrkW6jQeeCnqOKtQn67AyPB5K+CjCv7epdXnlGCd0upz0qNW/l5kAFsJFuNIhXi6A2uA5uH2U8BVUccVxjIUWAxkAZnAm0D/iGP61L0X+BVwS/j8FuDuFIrtGGAQkAOMTsH37iwgM3x+d4q9d63jnn+9LI+L8qGW8E8bA6x099XuXgQ8AUwsV2Yi8FD4fCpwuplZEmOsikTqk1bc/R1g92GKTAQe9sB0oK2ZdU1OdNWTQJ3Sirtvcfe54fP9wDKC5CBeWn1OCdZJGp7TgVXunkoremYCzc0skyDh3RxxPGWOIfjP9gF3LwH+C1wYZUCV3Hvjv+MfAs5PZkxlKorN3Ze5+/Io4imvkvheDz9bgOlAj6QHRqWx7YvbbAFEPjOJkvBP6w5siNveyKe/aA+VCf+y5QLtkxJd1SVSH4Avhj/TTDWznskJrc4kWud0M87MFpjZv81sSNTBJCrsrnUcMKPcobT9nA5TJ0jTz0mqbRLweNRBlHH3TcBvgPXAFiDX3V+PNqpDFgOfMbP2ZpZF8GtYKn7fdHb3LeHzrUDnKINJY9cQ/NqZMszsTjPbAFwORN5dUEm4ALwI9HH3Y4E3+LgFQFLHXIKfu4cDfwaejzacxJhZS+AZ4JvlWiHS1hHqlJafk1SPmTUBzgOejjqWMmH/5YlAX6Ab0MLMvhxtVAF3X0bQReF14FVgPlAaZUxH4kHfhchbTNONmd0GlACPRR1LPHe/zd17EsR1c9TxKAn/tE188n/mPcJ9FZYJf+5rA+xKSnRVd8T6uPsudz8Ybt4PjEpSbHUlkc8wrbj7PnfPC5+/AjQ2sw4Rh3VYZtaYIFl9zN2fraBI2n1OR6pTOn5OUiOfA+a6+7aoA4lzBrDG3Xe4ezHwLHBixDEd4u4PuPsodz8F2EMwtiLVbCvrGhf+uT3ieNKKmV0FfB64PPxPTCp6DPhi1EEoCf+0WcAAM+sbtnJMAqaVKzMNmBw+vwh4K4X/oh2xPuX64Z5H0Nc1nU0Drgxn3ziB4OfYLUc6KZWZWZeycQdmNobg326q/sePMNYHgGXu/rtKiqXV55RIndLtc5Iau5QU6ooSWg+cYGZZ4d/F00mhe7qZdQr/7EXQH/xf0UZUofjv+MnACxHGklbMbALwPeA8dz8QdTzxzGxA3OZE4MOoYimTGXUAqcbdS8zsZuA1glHvU9x9iZndAcx292kEX8SPmNlKgo7/k6KL+PASrM/Xzew8gp+OdhPMlpKyzOxxglkoOpjZRuDHQGMAd/8b8ApBX8OVwAHg6mgiTVwCdboIuMnMSoACYFIK/8cP4CTgCmCRmc0P9/0A6AVp+zklUqd0+5ykmsysBXAm8JWoY4nn7jPMbCpB16gSYB6ptVz3M2bWHigGvubue6MMppJ7713AU2Z2LbAO+FIKxbaboKtbR+BlM5vv7menUHy3Ak2BN8L2iOnufmOKxHaOmQ0CYgSfa9LjKk/L1ouIiIiIJJm6o4iIiIiIJJmScBERERGRJFMSLiIiIiKSZErCRURERESSTEm4iEgtMrMpZrbdzBYnUPYUM5trZiVmdlHc/t7h/vlmtsTMIh/FLyIitUuzo4iI1CIzOwXIAx5296FHKNsHaA18B5jm7lPD/U0I7s8HwxU6FwMnuvvmOg1eRESSRi3hIiK1yN3fIZjL9xAz62dmr5rZHDN718yODsuudfeFBPPWxl+jKG4V26boXi0iUu/oxi5SQ2Y21sw+MLMCM9tjZj+KOiZJOfcB/+fuowhave890glm1tPMFgIbgLvVCi4NhZk9aGYvRR0HgJllm9k2M+tXS9d72sy+XRvXkvSnJFykBszsDOBlglVUhwO/Au4ws5GRBiYpI+xOciLwdLjS5t+Brkc6z903uPuxQH9gspl1rtNARaQiPwBecfdVtXS9O4DbzKxNLV1P0piScJFqCvvt/gP4rrvf7+4fufsvga3AeDM718z+Em2UkgIaAXvdfUTc45hETw5bwBcDn6mzCEXkU8wsC7iOoJGlVrj7ImA18OXauqakLyXhItV3KtAWeLTc/mLgIHAsMD+5IUmqcfd9wBozuxjAAsMPd46Z9TCz5uHzbOBkYHmdByuSYsysqZn9IewSUmhm083s5HJlWpjZw2aWF5a71cxeMrMHa/jy5wAO/K+CuGrSDXEacGkNY5N6QEm4SPV9Fljk7sVlO8ysE9AdmEuQhB8dDsZbWjYYT+o3M3sc+AAYZGYbzexa4HLgWjNbACwBJoZljzezjcDFwN/NbEl4mWOAGWH5/wK/CVvQRBqaXwGXANcAxwGLgFfNLL5L128JGkUuILgvD6d2fjn6DDDHy00jVwvdEGcCY8r+oy0Nl6YoFKkmM3sVyHb3sXH7fkKQcA0k6EJwn7v/0cxuAMa4+3WRBCsikibCFuwOBMn3HuA6d384PJYBfAQ87u4/DMdc7AaudPcnwjItgI3AC+5+VbjvOWA88B93j5+T//MESXwjggHQ98cdex7IdffJcfuaEPwqdYe7/zNu/xbg1+7+OzM7F/icu99cSf2OBRYA/Wuxr7mkIbWEi1TfcQStndea2UAz+y7wfeBqoAmQBfw5LDuf4EtFREQS0w9oTFx3EHcvJfilaXC5MjPjyuQTNILE+yNwZfwOM8sEfkfQen4c8F0zax9XpDlQWO46R+qGCEfuilgQd31pwJSEi1SDmXUDOgGXAf9H8BPpZcBEd3+P4AtimbuXzf88ElgYRawiIvVQlX7Gd/ccYH+53WOAJe6+yd3zgH8DZ8Ud3wlklzvnSN0Q4chdEduFf+6oSh2k/lESLlI9IwhmvHglnO2iqbsf5+6vh8ePBfqZWePwBn0dH7eKi4jIka0CioCTynaE3VHGAUvjyhQDx8eVyQIOu1ptqBuwKW57E0EyXWYeH7e4lzmOYAGteF8lmPFkerh9LLAhXBfgDwRrA8QbCmxy920JxCj1WGbUAYikqeM4fMv2scBLwCwgA/h/7q5WDxGRBLl7vpn9FbjbzHYCa4BvAZ0JF7xy9zwzmxJXZgvwQ4JGxpoOenstvG57d98V7jsOaBoOuH6XYJD194Ez3d3NrCmf7op4Trnrfia8tjRwSsJFquewSbi7a0U0EZGa+3745z8J+mLPAya4+5a4Mt8BWhBM/ZcH/J4gUS/fn7u8zXyy5bs7n+xbvsjMZgKTgHviuiGeC/yC4D8CS/m4GyIcoSuimTUjmMXl7CNVXOo/zY4iIiIi9UbYGr2OYLaS38btHw/cXDY7SjgwcxnBrCm5wBzgxLhWb8xsAsGgzsEEifNj7l6+n3j8a08mWGVzKEF/8lcIZkrZER7/GkHSflZl15CGQy3hIiIikrbM7DiCufVnAq0IWs9bAU/GlXmTYE7vFmVz87v7B2b2beBtgu4rv4pPwAHc/VUzuwfowZG7IcKRuyIWEwzmF1FLuIiIiKSvMAn/BzAIKCHoh/0dd59Ty68zFdji7kqipVYoCRcRERERSTJNUSgiIiIikmRKwkVEREREkkxJuIiIiIhIkikJFxERERFJMiXhIiIiIiJJpiRcRERERCTJlISLiIiIiCTZ/wcspuipsPh29wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Set rho in [rho_ppoly_tab[0]/10,rho_ppoly_tab[2]*10]\n", "rho_for_plot = np.linspace(rho_ppoly_tab[0]/10,rho_ppoly_tab[NEOS-2]*10,1000)\n", "\n", "# Set P_{cold}(rho)\n", "eps_for_plot = np.linspace(0,0,1000)\n", "for i in range(len(eps_for_plot)):\n", " eps_for_plot[i] = eps_cold(rho_for_plot[i],rho_ppoly_tab,Gamma_ppoly_tab,K_ppoly_tab,NEOS)\n", "\n", "# Plot P_{cold}(rho) x rho\n", "import matplotlib.pyplot as plt\n", "f = plt.figure(figsize=(12,4))\n", "ax1 = f.add_subplot(121)\n", "ax1.set_title(r\"Plot #1: $\\varepsilon_{\\rm cold}\\times\\rho_{b}$\",fontsize='16')\n", "ax1.set_xlabel(r\"$\\rho_{b}$\",fontsize='14')\n", "ax1.set_ylabel(r\"$\\varepsilon_{\\rm cold}$\",fontsize='14')\n", "ax1.grid()\n", "ax1.plot(rho_for_plot,eps_for_plot)\n", "\n", "ax2 = f.add_subplot(122)\n", "ax2.set_title(r\"Plot #2: $\\log_{10}\\left(\\varepsilon_{\\rm cold}\\right)\\times\\log_{10}\\left(\\rho_{b}\\right)$\",fontsize='16')\n", "ax2.set_xlabel(r\"$\\log_{10}\\left(\\rho_{b}\\right)$\",fontsize='14')\n", "ax2.set_ylabel(r\"$\\log_{10}\\left(\\varepsilon_{\\rm cold}\\right)$\",fontsize='14')\n", "ax2.grid()\n", "ax2.plot(np.log10(rho_for_plot),np.log10(eps_for_plot))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 5: The piecewise polytrope EOS C code \\[Back to [top](#toc)\\]\n", "$$\\label{ppeos__c_code}$$\n", "\n", "Now we will begin implementing a C code to handle all the topics we discussed here.\n", "\n", "\n", "\n", "## Step 5.a: Preliminary treatment of the input \\[Back to [top](#toc)\\]\n", "$$\\label{ppeos__c_code__prelim}$$\n", "\n", "Here we assume that we get as input the following values:\n", "\n", "1. neos: the number of polytropic EOSs we will be handling (equivalent to $N+1$ in the discussion above)\n", "1. rho_ppoly_tab: an array containing the values $\\left\\{\\rho_{0},\\rho_{1},\\ldots,\\rho_{{\\rm neos}-1}\\right\\}$. Notice that for both neos=1 and neos=2, the size of this array is $1$.\n", "1. K_ppoly_tab: an array of size neos containing only the value of $K_{0}$ in the entry K_ppoly_tab[0].\n", "1. Gamma_ppoly_tab: an array of size neos containing the values $\\left\\{\\Gamma_{0},\\Gamma_{1},\\ldots,\\Gamma_{N}\\right\\}$." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:06.582124Z", "iopub.status.busy": "2021-05-17T23:43:06.581486Z", "iopub.status.idle": "2021-05-17T23:43:06.583600Z", "shell.execute_reply": "2021-05-17T23:43:06.584019Z" } }, "outputs": [], "source": [ "import os,sys\n", "IGM_src_dir = os.path.join(\"..\",\"src\")\n", "Piecewise_Polytrope_EOS__C = os.path.join(IGM_src_dir,\"IllinoisGRMHD_EoS_lowlevel_functs.C\")" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:06.588988Z", "iopub.status.busy": "2021-05-17T23:43:06.588240Z", "iopub.status.idle": "2021-05-17T23:43:06.591371Z", "shell.execute_reply": "2021-05-17T23:43:06.591812Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing ../src/IllinoisGRMHD_EoS_lowlevel_functs.C\n" ] } ], "source": [ "%%writefile $Piecewise_Polytrope_EOS__C\n", "#ifndef ILLINOISGRMHD_EOS_FUNCTS_C_\n", "#define ILLINOISGRMHD_EOS_FUNCTS_C_\n", "/* Function : setup_K_ppoly_tab__and__eps_integ_consts()\n", " * Authors : Leo Werneck\n", " * Description : For a given set of EOS inputs, determine\n", " * values of K_ppoly_tab that will result in a\n", " * everywhere continuous P_cold function.\n", " * Dependencies: none\n", " *\n", " * Inputs : eos - a struct containing the following\n", " * relevant quantities:\n", " * : neos - number of polytropic EOSs used\n", " * : rho_ppoly_tab - array of rho values that determine\n", " * the polytropic EOS to be used.\n", " * : Gamma_ppoly_tab - array of Gamma_cold values to be\n", " * used in each polytropic EOS.\n", " * : K_ppoly_tab - array of K_ppoly_tab values to be used\n", " * in each polytropic EOS. Only K_ppoly_tab[0]\n", " * is known prior to the function call\n", " * : eps_integ_const - array of C_{j} values, which are the\n", " * integration constants that arrise when\n", " * determining eps_{cold} for a piecewise\n", " * polytropic EOS. This array should be\n", " * uninitialized or contain absurd values\n", " * prior to this function call.\n", " *\n", " * Outputs : K_ppoly_tab - fully populated array of K_ppoly_tab\n", " * to be used in each polytropic EOS.\n", " * : eps_integ_const - fully populated array of C_{j}'s,\n", " * used to compute eps_cold for\n", " * a piecewise polytropic EOS.\n", " */\n", "static void setup_K_ppoly_tab__and__eps_integ_consts(eos_struct &eos){" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 5.a.i: Determining $\\left\\{K_{1},K_{2},\\ldots,K_{\\rm neos}\\right\\}$ \\[Back to [top](#toc)\\]\n", "$$\\label{ppeos__c_code__prelim__computing_ktab}$$\n", "\n", "We start by computing the values $\\left\\{K_{1},K_{2},\\ldots,K_{\\rm neos}\\right\\}$ from the input values. Remember that the values of $K_{j}$ are obtained by demanding that the $P_{\\rm cold}$ function be everywhere continuous, resulting in the relation (see [the discussion above for the details](#introduction))\n", "\n", "$$\n", "\\boxed{K_{j} = K_{j-1}\\rho_{j-1}^{\\Gamma_{j-1}-\\Gamma_{j}}\\ ,\\ j\\in\\left[1,N\\right]}\\ .\n", "$$\n", "\n", "**Implementation note**: Because the case neos=1 is handled first, the attentive reader will notice that we already set $C_{0}=0$ below, *before* discussing its implementation. For more on that, please refer to the [next subsection, where we discuss determining $\\left\\{C_{0},C_{1},C_{2},\\ldots,C_{\\rm neos}\\right\\}$](#ppeos__c_code__prelim__computing_eps_integ_consts)." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:06.597151Z", "iopub.status.busy": "2021-05-17T23:43:06.596332Z", "iopub.status.idle": "2021-05-17T23:43:06.599453Z", "shell.execute_reply": "2021-05-17T23:43:06.599874Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/IllinoisGRMHD_EoS_lowlevel_functs.C\n" ] } ], "source": [ "%%writefile -a $Piecewise_Polytrope_EOS__C\n", "\n", " /* When neos = 1, we will only need the value K_ppoly_tab[0] and eps_integ_const[0].\n", " * Since our only polytropic EOS is given by\n", " * -----------------------------------\n", " * | P_{0} = K_{0} * rho ^ (Gamma_{0}) | ,\n", " * -----------------------------------\n", " * and, therefore,\n", " * ---------------------------------------------------------------\n", " * | eps_{0} = K_{0} * rho ^ (Gamma_{0}-1) / (Gamma_{0}-1) + C_{0} | ,\n", " * ---------------------------------------------------------------\n", " * we only need to set up K_{0} := K_ppoly_tab[0] and C_{0} := eps_integ_const[0].\n", " * K_{0} is a user input, so we need to do nothing. C_{0}, on the other hand,\n", " * is fixed by demanding that eps(rho) -> 0 as rho -> 0. Thus, C_{0} = 0.\n", " */\n", " eos.eps_integ_const[0] = 0.0;\n", " if(eos.neos==1) return;\n", "\n", " /********************\n", " * Setting up K_{j} *\n", " ********************/\n", " /* When neos > 1, we have the following structure\n", " *\n", " * / K_{0} * rho^(Gamma_{0}) , rho < rho_{0}\n", " * | K_{1} * rho^(Gamma_{1}) , rho_{0} <= rho < rho_{1}\n", " * | ...\n", " * P(rho) = < K_{j} * rho^(Gamma_{j}) , rho_{j-1} <= rho < rho_{j}\n", " * | ...\n", " * | K_{neos-2} * rho^(Gamma_{neos-2}) , rho_{neos-3} <= rho < rho_{neos-2}\n", " * \\ K_{neos-1} * rho^(Gamma_{neos-1}) , rho >= rho_{neos-2}\n", " *\n", " * Imposing that P(rho) be everywhere continuous, we have\n", " * -------------------------------------------------\n", " * | K_{j} = K_{j-1} * rho^(Gamma_{j-1} - Gamma_{j}) |\n", " * -------------------------------------------------\n", " */\n", " for(int j=1; j\n", "\n", "### Step 5.a.ii: Determining $\\left\\{C_{0},C_{1},C_{2},\\ldots,C_{\\rm neos}\\right\\}$ \\[Back to [top](#toc)\\]\n", "$$\\label{ppeos__c_code__prelim__computing_eps_integ_consts}$$\n", "\n", "We now focus our attention on the computation of $\\left\\{C_{0},C_{1},C_{2},\\ldots,C_{\\rm neos}\\right\\}$, which are the integration constants that emerge when computing $\\epsilon_{\\rm cold}$ for a piecewise polytrope EOS (see [the discussion above for the details](#eps_cold)). Here we set them so that $\\epsilon_{\\rm cold}$ is also everywhere continuous, i.e.\n", "\n", "$$\n", "\\boxed{\n", "C_{j} =\n", "\\left\\{\n", "\\begin{matrix}\n", "0 & , & {\\rm if}\\ j=0\\\\\n", "C_{j-1} + \\frac{K_{j-1}\\rho_{j-1}^{\\Gamma_{j-1}-1}}{\\Gamma_{j-1}-1} - \\frac{K_{j}\\rho_{j-1}^{\\Gamma_{j}-1}}{\\Gamma_{j}-1} & , & {\\rm otherwise}\n", "\\end{matrix}\n", "\\right.\n", "}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:06.604403Z", "iopub.status.busy": "2021-05-17T23:43:06.603643Z", "iopub.status.idle": "2021-05-17T23:43:06.606487Z", "shell.execute_reply": "2021-05-17T23:43:06.606906Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/IllinoisGRMHD_EoS_lowlevel_functs.C\n" ] } ], "source": [ "%%writefile -a $Piecewise_Polytrope_EOS__C\n", "\n", " /********************\n", " * Setting up C_{j} *\n", " ********************/\n", " /* When neos > 1, we have the following structure (let neos->N):\n", " *\n", " * / K_{0}*rho^(Gamma_{0}-1)/(Gamma_{0}-1) + C_{0}, rho < rho_{0}\n", " * | K_{1}*rho^(Gamma_{1}-1)/(Gamma_{1}-1) + C_{1}, rho_{0} <= rho < rho_{1}\n", " * | ...\n", " * eps(rho) = < K_{j}*rho^(Gamma_{j}-1)/(Gamma_{j}-1) + C_{j}, rho_{j-1} <= rho < rho_{j}\n", " * | ...\n", " * | K_{N-2}*rho^(Gamma_{N-2}-1)/(Gamma_{N-2}-1) + C_{N-2}, rho_{N-3} <= rho < rho_{N-2}\n", " * \\ K_{N-1}*rho^(Gamma_{N-1}-1)/(Gamma_{N-1}-1) + C_{N-1}, rho >= rho_{N-2}\n", " *\n", " * Imposing that eps_{cold}(rho) be everywhere continuous, we have\n", " * ---------------------------------------------------------------\n", " * | C_{j} = C_{j-1} |\n", " * | + ( K_{j-1}*rho_{j-1}^(Gamma_{j-1}-1) )/(Gamma_{j-1}-1) |\n", " * | - ( K_{j+0}*rho_{j-1}^(Gamma_{j+0}-1) )/(Gamma_{j+0}-1) |\n", " * ---------------------------------------------------------------\n", " */\n", " for(int j=1; j\n", "\n", "## Step 5.b: Setting up the `eos_struct` \\[Back to [top](#toc)\\]\n", "$$\\label{ppeos__c_code__eos_struct_setup}$$\n", "\n", "We will now set up the `eos_struct` from the input. The `eos_struct` declaration can be found inside the `IllinoisGRMHD_headers.h` source file, but we will repeat it here for the sake of the reader:\n", "\n", "```c\n", "#define MAX_EOS_PARAMS 10\n", "struct eos_struct {\n", " int neos;\n", " CCTK_REAL rho_ppoly_tab[MAX_EOS_PARAMS-1];\n", " CCTK_REAL eps_integ_const[MAX_EOS_PARAMS],K_ppoly_tab[MAX_EOS_PARAMS],Gamma_ppoly_tab[MAX_EOS_PARAMS];\n", " CCTK_REAL Gamma_th;\n", "};\n", "```" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:06.611781Z", "iopub.status.busy": "2021-05-17T23:43:06.610859Z", "iopub.status.idle": "2021-05-17T23:43:06.614259Z", "shell.execute_reply": "2021-05-17T23:43:06.614963Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/IllinoisGRMHD_EoS_lowlevel_functs.C\n" ] } ], "source": [ "%%writefile -a $Piecewise_Polytrope_EOS__C\n", "\n", "/* Function : initialize_EOS_struct_from_input()\n", " * Authors : Leo Werneck\n", " * Description : Initialize the eos struct from user\n", " * input\n", " * Dependencies: setup_K_ppoly_tab__and__eps_integ_consts()\n", " * : cctk_parameters.h (FIXME)\n", " *\n", " * Inputs : eos - a struct containing the following\n", " * relevant quantities:\n", " * : neos - number of polytropic EOSs used\n", " * : rho_ppoly_tab - array of rho values that determine\n", " * the polytropic EOS to be used.\n", " * : Gamma_ppoly_tab - array of Gamma_cold values to be\n", " * used in each polytropic EOS.\n", " * : K_ppoly_tab - array of K_ppoly_tab values to be used\n", " * in each polytropic EOS.\n", " * : eps_integ_const - array of C_{j} values, which are the\n", " * integration constants that arrise when\n", " * determining eps_{cold} for a piecewise\n", " * polytropic EOS.\n", " *\n", " * Outputs : eos - fully initialized EOS struct\n", " */\n", "static void initialize_EOS_struct_from_input(eos_struct &eos){\n", " /* We start by setting up the eos_struct\n", " * with the inputs given by the user at\n", " * the start of the simulation. Keep in\n", " * mind that these parameters are found\n", " * in the \"cctk_parameters.h\" header file.\n", " * ^^^^^^^FIXME^^^^^^^\n", " */\n", "\n", " // Initialize: neos, {rho_{j}}, {K_{0}}, and {Gamma_{j}}\n", "#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", " DECLARE_CCTK_PARAMETERS;\n", "#endif\n", "\n", " eos.neos = neos;\n", " eos.K_ppoly_tab[0] = K_ppoly_tab0;\n", " for(int j=0; j<=neos-2; j++) eos.rho_ppoly_tab[j] = rho_ppoly_tab_in[j];\n", " for(int j=0; j<=neos-1; j++) eos.Gamma_ppoly_tab[j] = Gamma_ppoly_tab_in[j];\n", "\n", " // Initialize {K_{j}}, j>=1, and {eps_integ_const_{j}}\n", " setup_K_ppoly_tab__and__eps_integ_consts(eos);\n", "}\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 5.c: The `find_polytropic_K_and_Gamma_index()` function \\[Back to [top](#toc)\\]\n", "$$\\label{ppeos__c_code__find_polytropic_k_and_gamma_index}$$" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:06.619676Z", "iopub.status.busy": "2021-05-17T23:43:06.618884Z", "iopub.status.idle": "2021-05-17T23:43:06.621776Z", "shell.execute_reply": "2021-05-17T23:43:06.622187Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/IllinoisGRMHD_EoS_lowlevel_functs.C\n" ] } ], "source": [ "%%writefile -a $Piecewise_Polytrope_EOS__C\n", "\n", "/* Function : find_polytropic_K_and_Gamma_index()\n", " * Authors : Leo Werneck & Zach Etienne\n", " * Description : For a given value of rho, find the\n", " * appropriate values of Gamma_ppoly_tab\n", " * and K_ppoly_tab by determining the appropriate\n", " * index\n", " * Dependencies: initialize_EOS_struct_from_input()\n", " * : cctk_parameters.h (FIXME)\n", " *\n", " * Inputs : eos - a struct containing the following\n", " * relevant quantities:\n", " * : neos - number of polytropic EOSs used\n", " * : rho_ppoly_tab - array of rho values that determine\n", " *\n", " * Outputs : index - the appropriate index for the K_ppoly_tab\n", " * and Gamma_ppoly_tab array\n", " */\n", "static inline int find_polytropic_K_and_Gamma_index(eos_struct eos, CCTK_REAL rho_in) {\n", "\n", " /* We want to find the appropriate polytropic EOS for the\n", " * input value rho_in. Remember that:\n", " *\n", " * if rho < rho_{0}: P_{0} , index: 0\n", " * if rho >= rho_{0} but < rho_{1}: P_{1} , index: 1\n", " * if rho >= rho_{1} but < rho_{2}: P_{2} , index: 2\n", " * ...\n", " * if rho >= rho_{j-1} but < rho_{j}: P_{j} , index: j\n", " *\n", " * Then, a simple way of determining the index is through\n", " * the formula:\n", " * ---------------------------------------------------------------------------\n", " * | index = (rho >= rho_{0}) + (rho >= rho_{1}) + ... + (rho >= rho_{neos-2}) |\n", " * ---------------------------------------------------------------------------\n", " */\n", " if(eos.neos == 1) return 0;\n", "\n", " int polytropic_index = 0;\n", " for(int j=0; j<=eos.neos-2; j++) polytropic_index += (rho_in >= eos.rho_ppoly_tab[j]);\n", "\n", " return polytropic_index;\n", "}\n", "\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:06.626947Z", "iopub.status.busy": "2021-05-17T23:43:06.626080Z", "iopub.status.idle": "2021-05-17T23:43:06.630053Z", "shell.execute_reply": "2021-05-17T23:43:06.629606Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/IllinoisGRMHD_EoS_lowlevel_functs.C\n" ] } ], "source": [ "%%writefile -a $Piecewise_Polytrope_EOS__C\n", "\n", "/* Function : compute_P_cold__eps_cold()\n", " * Authors : Leo Werneck\n", " * Description : Computes P_cold and eps_cold.\n", " * Dependencies: initialize_EOS_struct_from_input()\n", " * : find_polytropic_K_and_Gamma_index()\n", " *\n", " * Inputs : P_cold - cold pressure\n", " * : eps_cold - cold specific internal energy\n", " * : eos - a struct containing the following\n", " * relevant quantities:\n", " * : neos - number of polytropic EOSs used\n", " * : rho_ppoly_tab - array of rho values that determine\n", " * the polytropic EOS to be used.\n", " * : Gamma_ppoly_tab - array of Gamma_cold values to be\n", " * used in each polytropic EOS.\n", " * : K_ppoly_tab - array of K_ppoly_tab values to be used\n", " * in each polytropic EOS.\n", " * : eps_integ_const - array of C_{j} values, which are the\n", " * integration constants that arrise when\n", " * determining eps_{cold} for a piecewise\n", " * polytropic EOS.\n", " *\n", " * Outputs : P_cold - cold pressure (supports SPEOS and PPEOS)\n", " * : eps_cold - cold specific internal energy (supports SPEOS and PPEOS)\n", " * : polytropic_index - polytropic index used for P_cold and eps_cold\n", " *\n", " * SPEOS: Single-Polytrope Equation of State\n", " * PPEOS: Piecewise Polytrope Equation of State\n", " */\n", "static inline void compute_P_cold__eps_cold(eos_struct eos, CCTK_REAL rho_in,\n", " CCTK_REAL &P_cold,CCTK_REAL &eps_cold) {\n", " // This code handles equations of state of the form defined\n", " // in Eqs 13-16 in http://arxiv.org/pdf/0802.0200.pdf\n", " if(rho_in==0) {\n", " P_cold = 0.0;\n", " eps_cold = 0.0;\n", " return;\n", " }\n", "\n", " /* --------------------------------------------------\n", " * | Single and Piecewise Polytropic EOS modification |\n", " * --------------------------------------------------\n", " *\n", " * We now begin our modifications to this function so that\n", " * it supports both single and piecewise polytropic equations\n", " * of state.\n", " *\n", " * The modifications below currently assume that the user\n", " * has called the recently added function\n", " *\n", " * - initialize_EOS_struct_from_input()\n", " *\n", " * *before* this function is called. We can add some feature\n", " * to check this automatically as well, but we'll keep that as\n", " * a TODO/FIXME for now.\n", " */\n", "\n", " /* First, we compute the pressure, which in the case of a\n", " * piecewise polytropic EOS is given by\n", " *\n", " * / P_{1} / K_{1} * rho^(Gamma_{1}) , rho_{0} <= rho < rho_{1}\n", " * | ... | ...\n", " * P(rho) = < P_{j} = < K_{j} * rho^(Gamma_{j}) , rho_{j-1} <= rho < rho_{j}\n", " * | ... | ...\n", " * \\ P_{neos-2} \\ K_{neos-2} * rho^(Gamma_{neos-2}), rho_{neos-3} <= rho < rho_{neos-2}\n", " *\n", " * The index j is determined by the find_polytropic_K_and_Gamma_index() function.\n", " */\n", " // Set up useful auxiliary variables\n", " int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_in);\n", " CCTK_REAL K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];\n", " CCTK_REAL Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];\n", " CCTK_REAL eps_integ_const = eos.eps_integ_const[polytropic_index];\n", "\n", " // Then compute P_{cold}\n", " P_cold = K_ppoly_tab*pow(rho_in,Gamma_ppoly_tab);\n", "\n", " /* Then we compute the cold component of the specific internal energy,\n", " * which in the case of a piecewise polytropic EOS is given by (neos -> N)\n", " *\n", " * / P_{1}/(rho*(Gamma_{1}-1)) + C_{1} , rho_{0} <= rho < rho_{1}\n", " * | ...\n", " * eps(rho) = < P_{j}/(rho*(Gamma_{j}-1)) + C_{j} , rho_{j-1} <= rho < rho_{j}\n", " * | ...\n", " * \\ P_{N-2}/(rho*(Gamma_{N-2}-1)) + C_{N-2}, rho_{N-3} <= rho < rho_{N-2}\n", " */\n", " eps_cold = P_cold/(rho_in*(Gamma_ppoly_tab-1.0)) + eps_integ_const;\n", "\n", "}\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 5.d: The new `compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold()` function \\[Back to [top](#toc)\\]\n", "$$\\label{ppeos__c_code__compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold}$$\n", "\n", "We now write down the updated version of the `compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold()` function, which is found in the `inlined_functions.C` source file of `IllinoisGRMHD`, and documented in the [`IllinoisGRMHD`'s inlined functions NRPy+ tutorial module](Tutorial-IllinoisGRMHD__inlined_functions.ipynb)." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:06.634785Z", "iopub.status.busy": "2021-05-17T23:43:06.634017Z", "iopub.status.idle": "2021-05-17T23:43:06.636824Z", "shell.execute_reply": "2021-05-17T23:43:06.637230Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/IllinoisGRMHD_EoS_lowlevel_functs.C\n" ] } ], "source": [ "%%writefile -a $Piecewise_Polytrope_EOS__C\n", "\n", "/* Function : compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold()\n", " * Authors : Leo Werneck & Zach Etienne\n", " * Description : Compute basic quantities related to\n", " * : the EOS, namely: P_cold, eps_cold,\n", " * : dPcold/drho, eps_th, h, and Gamma_cold\n", " * Dependencies: initialize_EOS_struct_from_input()\n", " *\n", " * Inputs : U - array containing primitives {rho,P,v^{i},B^i}\n", " * : P_cold - cold pressure\n", " * : eps_cold - cold specific internal energy\n", " * : dPcold_drho - derivative of P_cold\n", " * : eps_th - thermal specific internal energy\n", " * : h - enthalpy\n", " * : Gamma_cold - cold polytropic Gamma\n", " * : eos - a struct containing the following\n", " * relevant quantities:\n", " * : neos - number of polytropic EOSs used\n", " * : rho_ppoly_tab - array of rho values that determine\n", " * the polytropic EOS to be used.\n", " * : Gamma_ppoly_tab - array of Gamma_cold values to be\n", " * used in each polytropic EOS.\n", " * : K_ppoly_tab - array of K_ppoly_tab values to be used\n", " * in each polytropic EOS.\n", " * : eps_integ_const - array of C_{j} values, which are the\n", " * integration constants that arrise when\n", " * determining eps_{cold} for a piecewise\n", " * polytropic EOS.\n", " *\n", " * Outputs : P_cold - cold pressure (supports SPEOS and PPEOS)\n", " * : eps_cold - cold specific internal energy (supports SPEOS and PPEOS)\n", " * : dPcold_drho - derivative of P_cold (supports SPEOS and PPEOS)\n", " * : eps_th - thermal specific internal energy (supports SPEOS and PPEOS)\n", " * : h - enthalpy (supports SPEOS and PPEOS)\n", " * : Gamma_cold - cold polytropic Gamma (supports SPEOS and PPEOS)\n", " *\n", " * SPEOS: Single-Polytrope Equation of State\n", " * PPEOS: Piecewise Polytrope Equation of State\n", " */\n", "static inline void compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold(CCTK_REAL *U, eos_struct &eos, CCTK_REAL Gamma_th,\n", " CCTK_REAL &P_cold,CCTK_REAL &eps_cold,CCTK_REAL &dPcold_drho,CCTK_REAL &eps_th,CCTK_REAL &h,\n", " CCTK_REAL &Gamma_cold) {\n", " // This code handles equations of state of the form defined\n", " // in Eqs 13-16 in http://arxiv.org/pdf/0802.0200.pdf\n", "\n", " if(U[RHOB]==0) {\n", " P_cold = 0.0;\n", " eps_cold = 0.0;\n", " dPcold_drho = 0.0;\n", " eps_th = 0.0;\n", " h = 0.0;\n", " Gamma_cold = eos.Gamma_ppoly_tab[0];\n", " return;\n", " }\n", "\n", " int polytropic_index = find_polytropic_K_and_Gamma_index(eos, U[RHOB]);\n", " compute_P_cold__eps_cold(eos,U[RHOB], P_cold,eps_cold);\n", " CCTK_REAL Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];\n", "\n", " // Set auxiliary variable rho_b^{-1}\n", " CCTK_REAL U_RHOB_inv = 1.0/U[RHOB];\n", "\n", " // Next compute dP/drho = Gamma * P / rho\n", " dPcold_drho = Gamma_ppoly_tab*P_cold*U_RHOB_inv;\n", "\n", " // Then we compute eps_th, h, and set Gamma_cold = Gamma_ppoly_tab[j].\n", " eps_th = (U[PRESSURE] - P_cold)/(Gamma_th-1.0)*U_RHOB_inv;\n", " h = 1.0 + eps_cold + eps_th + U[PRESSURE]*U_RHOB_inv;\n", " Gamma_cold = Gamma_ppoly_tab;\n", "\n", "}\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 5.e: The `print_EOS_table()` function \\[Back to [top](#toc)\\]\n", "$$\\label{ppeos__c_code__print_eos_table}$$\n", "\n", "This function, which is called for diagnostic purposes and for user convenience, prints out the most relevant EOS table components to the user at the *beginning* of the run. Notice that the function call happens only once, in the [`driver_conserv_to_prims.C`](Tutorial-IllinoisGRMHD__the_conservative_to_primitive_algorithm.ipynb) file." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:06.641481Z", "iopub.status.busy": "2021-05-17T23:43:06.640728Z", "iopub.status.idle": "2021-05-17T23:43:06.643385Z", "shell.execute_reply": "2021-05-17T23:43:06.643809Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/IllinoisGRMHD_EoS_lowlevel_functs.C\n" ] } ], "source": [ "%%writefile -a $Piecewise_Polytrope_EOS__C\n", "\n", "/* Function : print_EOS_table()\n", " * Authors : Leo Werneck\n", " * Description : Prints out the EOS table, for diagnostic purposes\n", " *\n", " * Dependencies: initialize_EOS_struct_from_input()\n", " *\n", " * Inputs : eos - a struct containing the following\n", " * relevant quantities:\n", " * : neos - number of polytropic EOSs used\n", " * : rho_ppoly_tab - array of rho values that determine\n", " * the polytropic EOS to be used.\n", " * : Gamma_ppoly_tab - array of Gamma_cold values to be\n", " * used in each polytropic EOS.\n", " *\n", " * Outputs : CCTK_VInfo string with the EOS table used by IllinoisGRMHD\n", " */\n", "static inline void print_EOS_table( eos_struct eos ) {\n", "\n", " /* Start by printint a header t the table */\n", "#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", " CCTK_VInfo(CCTK_THORNSTRING,\"\\n\"\n", "#else\n", " printf(\"\\n\"\n", "#endif\n", "\".--------------------------------------------.\\n\"\n", "\"| Piecewise Polytropic EOS Parameters |\\n\"\n", "\".--------------------------------------------.\");\n", "\n", " /* Adjust the maximum index of rhob to\n", " * allow for single polytropes as well\n", " */\n", " int max_rho_index;\n", " if( eos.neos==1 ) {\n", " max_rho_index = 0;\n", " }\n", " else {\n", " max_rho_index = eos.neos-1;\n", " }\n", "\n", " /* Print out rho_pppoly_tab */\n", " for(int jj=0; jj\n", "\n", "# Step n: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n", "$$\\label{latex_pdf_output}$$\n", "\n", "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n", "[Tutorial-IllinoisGRMHD__piecewise_polytrope_tests.pdf](Tutorial-IllinoisGRMHD__piecewise_polytrope_tests.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means)." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "execution": { "iopub.execute_input": "2021-05-17T23:43:06.648564Z", "iopub.status.busy": "2021-05-17T23:43:06.647948Z", "iopub.status.idle": "2021-05-17T23:43:06.771062Z", "shell.execute_reply": "2021-05-17T23:43:06.771834Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "zsh:1: no matches found: Tut*.out\n" ] } ], "source": [ "import os\n", "nrpy_dir_path = os.path.join(\"..\",\"..\")\n", "latex_nrpy_style_path = os.path.join(nrpy_dir_path,\"latex_nrpy_style.tplx\")\n", "#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__Piecewise_Polytrope_EOS.ipynb\n", "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__Piecewise_Polytrope_EOS.tex\n", "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__Piecewise_Polytrope_EOS.tex\n", "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__Piecewise_Polytrope_EOS.tex\n", "!rm -f Tut*.out Tut*.aux Tut*.log" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.11" } }, "nbformat": 4, "nbformat_minor": 4 }