diff --git a/examples/pinball-demo.ipynb b/examples/pinball-demo.ipynb index 5ca9351..8377152 100644 --- a/examples/pinball-demo.ipynb +++ b/examples/pinball-demo.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "64a55749", "metadata": {}, "outputs": [], @@ -14,7 +14,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "id": "6ced7049", "metadata": {}, "outputs": [], @@ -22,29 +22,32 @@ "from pinballrt.sources import Star\n", "from pinballrt.grids import UniformCartesianGrid\n", "from pinballrt.model import Model\n", + "from pinballrt.gas import Gas\n", "\n", "import matplotlib.pyplot as plt\n", "import astropy.units as u\n", "import numpy as np" ] }, - { - "cell_type": "code", - "execution_count": 2, - "id": "0d7c5715", - "metadata": {}, - "outputs": [], - "source": [ - "star = Star()\n", - "star.set_blackbody_spectrum()" - ] - }, { "cell_type": "code", "execution_count": 3, "id": "a8890f61", "metadata": {}, "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Warp CUDA error: Could not open libcuda.so.\n", + "Warp CUDA error: Function cuDriverGetVersion_f: a suitable driver entry point was not found\n", + "Warp CUDA error 36 (in function cuda_driver_version, /builds/omniverse/warp/warp/native/warp.cu:1719)\n", + "/opt/conda/envs/warp/lib/python3.13/site-packages/astropy/units/quantity.py:658: RuntimeWarning: overflow encountered in expm1\n", + " result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n", + "/opt/conda/envs/warp/lib/python3.13/site-packages/astropy/units/quantity.py:658: RuntimeWarning: overflow encountered in multiply\n", + " result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n" + ] + }, { "name": "stdout", "output_type": "stream", @@ -65,13 +68,6 @@ "name": "stderr", "output_type": "stream", "text": [ - "Warp CUDA error: Could not open libcuda.so.\n", - "Warp CUDA error: Function cuDriverGetVersion_f: a suitable driver entry point was not found\n", - "Warp CUDA error 36 (in function cuda_driver_version, /builds/omniverse/warp/warp/native/warp.cu:1719)\n", - "/opt/conda/envs/warp/lib/python3.13/site-packages/astropy/units/quantity.py:658: RuntimeWarning: overflow encountered in expm1\n", - " result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n", - "/opt/conda/envs/warp/lib/python3.13/site-packages/astropy/units/quantity.py:658: RuntimeWarning: overflow encountered in multiply\n", - " result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n", "/opt/conda/envs/warp/lib/python3.13/site-packages/astropy/units/quantity.py:658: RuntimeWarning: overflow encountered in expm1\n", " result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n", "/opt/conda/envs/warp/lib/python3.13/site-packages/astropy/units/quantity.py:658: RuntimeWarning: overflow encountered in multiply\n", @@ -82,15 +78,33 @@ "source": [ "model = Model(grid=UniformCartesianGrid, grid_kwargs={\"ncells\":9, \"dx\":2.0*u.au})\n", "\n", - "density = np.ones(model.grid.shape)*1.0e-16 * u.g / u.cm**3\n", + "density = np.ones(model.grid.shape)*1.0e-14 * u.g / u.cm**3\n", + "density[4,4,4] = 0.0 * u.g / u.cm**3 # Create a low-density cavity in the center\n", "\n", - "model.add_density(density, \"yso.dst\")\n", - "model.add_star(star)" + "vx, vy, vz = np.meshgrid(0.5*(model.grid.grid.w1.numpy()[1:] + model.grid.grid.w1.numpy()[0:-1]), \n", + " 0.5*(model.grid.grid.w2.numpy()[1:] + model.grid.grid.w2.numpy()[0:-1]), \n", + " 0.5*(model.grid.grid.w3.numpy()[1:] + model.grid.grid.w3.numpy()[0:-1]), indexing='ij')\n", + "velocity = np.concatenate((vx[np.newaxis], vy[np.newaxis], vz[np.newaxis]), axis=0) * (-1.0 * u.km / u.s)\n", + "\n", + "model.set_physical_properties(density=density, dust=\"yso.dst\", amax=1.0*u.mm, p=3.5, gases=['co.dat'], \n", + " abundances=[1.0e-4], microturbulence=0.2 * u.km / u.s, velocity=velocity)" ] }, { "cell_type": "code", "execution_count": 4, + "id": "7bc813aa", + "metadata": {}, + "outputs": [], + "source": [ + "star = Star()\n", + "star.set_blackbody_spectrum()\n", + "model.add_star(star)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, "id": "47dca2cd", "metadata": {}, "outputs": [ @@ -99,16 +113,18 @@ "output_type": "stream", "text": [ "Iteration 0\n", - "Module pinballrt.sources f5dedf3 load on device 'cpu' took 4.44 ms (cached)\n", - "Module pinballrt.grids 237c235 load on device 'cpu' took 10662.48 ms (compiled)\n", - "Module pinballrt.utils 08b2efb load on device 'cpu' took 1.45 ms (cached)\n" + "Module pinballrt.sources 1808836 load on device 'cpu' took 1.43 ms (cached)\n", + "Module pinballrt.grids 53af283 load on device 'cpu' took 0.46 ms (cached)\n", + "Module pinballrt.utils 08b2efb load on device 'cpu' took 0.29 ms (cached)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ - "100%|██████████| 1000000/1000000 [00:04<00:00, 227489.85it/s]\n" + "100%|██████████| 1000000/1000000 [00:06<00:00, 148369.58it/s]\n", + "/workspaces/pinball-warp/pinballrt/grids.py:466: RuntimeWarning: invalid value encountered in divide\n", + " temperature = ((total_energy*u.L_sun).cgs.value / (4*const.sigma_sb.cgs.value*\\\n" ] }, { @@ -123,7 +139,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "100%|██████████| 1000000/1000000 [00:43<00:00, 23251.28it/s]\n" + "100%|██████████| 1000000/1000000 [00:15<00:00, 63438.46it/s]\n" ] }, { @@ -138,14 +154,14 @@ "name": "stderr", "output_type": "stream", "text": [ - "100%|██████████| 1000000/1000000 [00:43<00:00, 23147.92it/s]\n" + "100%|██████████| 1000000/1000000 [00:48<00:00, 20424.26it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "2 1.1616004 230.01974\n", + "2 1.5704488 3.8647785\n", "Iteration 3\n" ] }, @@ -153,14 +169,44 @@ "name": "stderr", "output_type": "stream", "text": [ - "100%|██████████| 1000000/1000000 [00:52<00:00, 19069.58it/s]" + "100%|██████████| 1000000/1000000 [01:11<00:00, 14001.56it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3 1.1686034 1.343868\n", + "Iteration 4\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 1000000/1000000 [01:21<00:00, 12278.75it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "4 1.0577108 1.1048422\n", + "Iteration 5\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 1000000/1000000 [01:23<00:00, 11998.20it/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "3 1.192686 1.026761\n" + "5 1.0173565 1.0396658\n" ] }, { @@ -177,13 +223,13 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "id": "fa5b80fa", "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAGdCAYAAAAotLvzAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMERJREFUeJzt3XtwVHWa//HPIZdOYJJognToMmCYZUQJKoLLGFHCAnEjFx12ZRRHmdGtwuUiMQrKoLPRlbQwO5gpsuJi8QNGCuG3NYLM7qgEXUA2wxoC8YIW6JiFqGSy7o9NCEIndJ/fH0KvLcSk+3TS36bfr6pv1fTp8+3zpHHy5Hm+52LZtm0LAAAYo0+sAwAAAKFIzgAAGIbkDACAYUjOAAAYhuQMAIBhSM4AABiG5AwAgGFIzgAAGCa5tw8YCAT0xRdfKCMjQ5Zl9fbhAQAO2LatEydOyOPxqE+fnqvvTp8+rfb2dsefk5qaqrS0tChE1Lt6PTl/8cUXysvL6+3DAgCiqLGxUZdffnmPfPbp06eVP/h7amr2O/6s3NxcNTQ0xF2C7vXknJGRIUm6OfkOJVspvX34i4qVlhrrELpkJff6f2LhswOxjuCiYXeciXUIXbLPOP+F3+MC5t5V+YzdobfPbA3+Lu8J7e3tamr2q6FusDIzIq/OW08ElD/qiNrb20nOXTnXyk62UkjODllWHCRnKw6Ss0jO0WJb5p/GYlvm/wEhy9zkfE5vLEtmZvRxlJzjWTz85gQAJCC/HZDfwd8p/jjuipGcAQBGCshWQJFnZydzY43kDAAwUkABR4tOzmbHVmI28wEAMBiVMwDASH7blt+OvDXtZG6skZwBAEZK5DVn2toAABiGyhkAYKSAbPkTtHImOQMAjERbGwAAGIPKGQBgJM7WBgDAMAE5u/N9/N6ChLY2AADGiSg5P//888rPz1daWppGjRqlt99+O9pxAQASnP/s2dpORrwKOzlv3rxZpaWlWrJkiQ4cOKCbb75ZJSUlOnr0aE/EBwBIUH7b+YhXYSfnFStW6IEHHtDf/M3f6KqrrlJlZaXy8vK0atWqnogPAJCgAlEY8Sqs5Nze3q66ujoVFxeHbC8uLlZNTc0F5/h8PrW2toYMAADQubCS85dffim/3y+32x2y3e12q6mp6YJzvF6vsrKygiMvLy/yaAEACSMgS34HIyAr1j9CxCI6IcyyQn9g27bP23bO4sWL1dLSEhyNjY2RHBIAkGACtvMRr8K6zrl///5KSko6r0pubm4+r5o+x+VyyeVyRR4hAAAJJqzKOTU1VaNGjVJ1dXXI9urqahUWFkY1MABAYnPS0j434lXYdwgrKyvTvffeq9GjR+vGG2/U6tWrdfToUT344IM9ER8AIEE5TbAJlZx//OMf67//+7/19NNP69ixYyooKNDvf/97DR48uCfiAwAg4UR0b+05c+Zozpw50Y4FAICggG0pYEde/TqZG2s8+AIAYKREbmvz4AsAAAxD5QwAMJJffeR3UEP6oxhLbyM5AwCMZDtcc7ZZcwYAILpYcwYAAMagcgYAGMlv95HfdrDmnCj31gYAoLcEZCngoMEbUPxmZ9raAAAYhsoZAGAkTggDAMAw59acnYxw7d69W1OnTpXH45FlWdq6det5+3z00UeaNm2asrKylJGRoR/+8Ic6evRo8H2fz6f58+erf//+6tevn6ZNm6bPPvssrDhIzgAAnHXy5Elde+21qqqquuD7f/zjHzV27FgNGzZMO3fu1Lvvvqsnn3xSaWlpwX1KS0u1ZcsWbdq0SXv27FFbW5umTJkiv7/7t0WhrQ0AMNLXJ4Q5ePBFBHNLSkpUUlLS6ftLlizRbbfdpuXLlwe3DRkyJPi/W1patGbNGr300kuaOHGiJGnDhg3Ky8vTjh07dOutt3Yrjtgl5z6WZJm7HtAnPa3rnWLN5Yp1BF2yUlJiHULXDP7vMEQgEOsIutbeEesIumT5fLEOoUt2e3usQ+iUZVtSL/0zBxzevvPc2dqtra0h210ul1wR/P4MBAL613/9Vy1atEi33nqrDhw4oPz8fC1evFh33HGHJKmurk4dHR0qLi4OzvN4PCooKFBNTU23kzNtbQDARS0vL09ZWVnB4fV6I/qc5uZmtbW16dlnn9Vf/uVfavv27frRj36k6dOna9euXZKkpqYmpaam6tJLLw2Z63a71dTU1O1j0dYGABjJ+U1Ivq6cGxsblZmZGdweSdUsfV05S9Ltt9+uhx9+WJJ03XXXqaamRi+88ILGjRvX6VzbtmWF0aWjcgYAGCmgPo6HJGVmZoaMSJNz//79lZycrKuvvjpk+1VXXRU8Wzs3N1ft7e06fvx4yD7Nzc1yu93dPhbJGQBgJL9tOR7RlJqaqhtuuEGHDh0K2X748GENHjxYkjRq1CilpKSouro6+P6xY8f0wQcfqLCwsNvHoq0NAMBZbW1t+uSTT4KvGxoaVF9fr+zsbA0aNEgLFy7Uj3/8Y91yyy0aP368Xn/9df3ud7/Tzp07JUlZWVl64IEH9MgjjygnJ0fZ2dl69NFHNWLEiODZ291BcgYAGMnv8GxtfwT31t63b5/Gjx8ffF1WViZJmjVrltatW6cf/ehHeuGFF+T1evXQQw/pyiuv1G9/+1uNHTs2OOe5555TcnKyZsyYoVOnTmnChAlat26dkpKSuh2HZdt2r94ZvLW1VVlZWRrvmqFky9zLbLiUKjq4lCqK4uBSKjsOLqUSl1I5csZu11un/69aWlpCTrKKpnN54v/sH6m+Gd1PaN/21Qm/7r/+QI/G2lNYcwYAwDC0tQEARopFW9sUJGcAgJECkqMzrs1fCOocbW0AAAxD5QwAMNI3byQS6fx4RXIGABjJ+e074zc5x2/kAABcpKicAQBGisXznE1BcgYAGIm2dhh2796tqVOnyuPxyLIsbd26tQfCAgAkunPXOTsZ8SrsyE+ePKlrr71WVVVVPREPAAAJL+y2dklJiUpKSnoiFgAAggK2pYCTm5BE+ZGRvanH15x9Pp9837jRfGtra08fEgBwEQg4bE3H83XOPR651+tVVlZWcOTl5fX0IQEAiGs9npwXL16slpaW4GhsbOzpQwIALgIBu4/jEa96vK3tcrnkioPnDgMAzOKXJb+Da5WdzI21+P2zAgCAi1TYlXNbW5s++eST4OuGhgbV19crOztbgwYNimpwAIDE5bQ1nVBt7X379mn8+PHB12VlZZKkWbNmad26dVELDACQ2Pxy1pr2Ry+UXhd2ci4qKpJt2z0RCwAAEPfWBgAYirY2AACGSeQHX5CcAQBGsh0+MtLmUioAABAtVM4AACPR1gYAwDCJ/FSq+P2zAgCAixSVMwDASH6Hj4x0MjfWSM4AACPR1gYAAMagcgYAGCmgPgo4qCGdzI01kjMAwEh+25LfQWvaydxYi98/KwAAuEjFrHK2kpNkWQYX7impsY6gS1Z6WqxD6JKd7op1CF2yUw3+7/AbrA7zH4BnnW6PdQhdsvvEQTXlN/jf2k7qtUMl8glh8fFbCQCQcGyHT6Wy4/gOYfEbOQDgouaX5XiEa/fu3Zo6dao8Ho8sy9LWrVs73Xf27NmyLEuVlZUh230+n+bPn6/+/furX79+mjZtmj777LOw4iA5AwBw1smTJ3XttdeqqqrqO/fbunWr/uM//kMej+e890pLS7VlyxZt2rRJe/bsUVtbm6ZMmSJ/GMsVtLUBAEYK2M7WjQN2+HNKSkpUUlLynft8/vnnmjdvnt544w1Nnjw55L2WlhatWbNGL730kiZOnChJ2rBhg/Ly8rRjxw7deuut3YqDyhkAYKTA2TVnJ0OSWltbQ4bP54s8pkBA9957rxYuXKjhw4ef935dXZ06OjpUXFwc3ObxeFRQUKCamppuH4fkDAC4qOXl5SkrKys4vF5vxJ+1bNkyJScn66GHHrrg+01NTUpNTdWll14ast3tdqupqanbx6GtDQAwUkCWAhGc1PXN+ZLU2NiozMzM4HaXK7JLPOvq6vTrX/9a+/fvl2WFF5dt22HNoXIGABjp3B3CnAxJyszMDBmRJue3335bzc3NGjRokJKTk5WcnKwjR47okUce0RVXXCFJys3NVXt7u44fPx4yt7m5WW63u9vHIjkDANAN9957r9577z3V19cHh8fj0cKFC/XGG29IkkaNGqWUlBRVV1cH5x07dkwffPCBCgsLu30s2toAACMFHN6EJJK5bW1t+uSTT4KvGxoaVF9fr+zsbA0aNEg5OTkh+6ekpCg3N1dXXnmlJCkrK0sPPPCAHnnkEeXk5Cg7O1uPPvqoRowYETx7uztIzgAAIwXk8PadEaxX79u3T+PHjw++LisrkyTNmjVL69at69ZnPPfcc0pOTtaMGTN06tQpTZgwQevWrVNSUvdvfUpyBgDgrKKiItl29y+Q/s///M/ztqWlpWnlypVauXJlxHGQnAEARrIdnq1tO5gbayRnAICReCoVAACGicUJYaaI38gBALhIhZWcvV6vbrjhBmVkZGjAgAG64447dOjQoZ6KDQCQwM61tZ2MeBVWct61a5fmzp2rvXv3qrq6WmfOnFFxcbFOnjzZU/EBABLUudt3OhnxKqw159dffz3k9dq1azVgwADV1dXplltuiWpgAAAkKkcnhLW0tEiSsrOzO93H5/OFPJ6rtbXVySEBAAkikc/WjviEMNu2VVZWprFjx6qgoKDT/bxeb8ijuvLy8iI9JAAggbDmHIF58+bpvffe08svv/yd+y1evFgtLS3B0djYGOkhAQBICBG1tefPn69t27Zp9+7duvzyy79zX5fLFfHjuQAAiSuR29phJWfbtjV//nxt2bJFO3fuVH5+fk/FBQBIcCTnbpo7d642btyoV199VRkZGWpqapL09SOy0tPTeyRAAAASTVhrzqtWrVJLS4uKioo0cODA4Ni8eXNPxQcASFC2nF3r3P1nS5kn7LY2AAC9gbY2AACGSeTkzIMvAAAwDJUzAMBIiVw5k5wBAEZK5ORMWxsAAMNQOQMAjGTblmwH1a+TubFGcgYAGMnpM5nj+XnOtLUBADAMlTMAwEiJfEIYyRkAYKREXnOmrQ0AgGGonAEARqKtDQCAYRK5rR2z5GylJMuyzP3bwEpNiXUIXbL7mf8MbZ/7e7EOoUunBpj/by1J6c0dsQ6hS6nNJ2MdQpf6dJyJdQhdS/HFOoJOWXag145lO6yc4zk5s+YMAIBhzC1dAQAJzZZk287mxyuSMwDASAFZsrhDGAAAMAGVMwDASIl8tjaVMwDASOeuc3YywrV7925NnTpVHo9HlmVp69atwfc6Ojr02GOPacSIEerXr588Ho/uu+8+ffHFFyGf4fP5NH/+fPXv31/9+vXTtGnT9Nlnn4UVB8kZAICzTp48qWuvvVZVVVXnvffVV19p//79evLJJ7V//3698sorOnz4sKZNmxayX2lpqbZs2aJNmzZpz549amtr05QpU+T3+7sdB21tAICRbNvh2doRzC0pKVFJSckF38vKylJ1dXXItpUrV+rP//zPdfToUQ0aNEgtLS1as2aNXnrpJU2cOFGStGHDBuXl5WnHjh269dZbuxUHlTMAwEjn1pydjJ7W0tIiy7J0ySWXSJLq6urU0dGh4uLi4D4ej0cFBQWqqanp9udSOQMALmqtra0hr10ul1wul+PPPX36tB5//HHNnDlTmZmZkqSmpialpqbq0ksvDdnX7Xarqamp259N5QwAMFK0Kue8vDxlZWUFh9frdRxbR0eH7rrrLgUCAT3//PPd+FlsWVb3K3kqZwCAkQK2JSsKT6VqbGwMVraSHFfNHR0dmjFjhhoaGvTWW2+FfHZubq7a29t1/PjxkOq5ublZhYWF3T4GlTMAwEjnTghzMiQpMzMzZDhJzucS88cff6wdO3YoJycn5P1Ro0YpJSUl5MSxY8eO6YMPPggrOVM5AwBwVltbmz755JPg64aGBtXX1ys7O1sej0d//dd/rf379+tf/uVf5Pf7g+vI2dnZSk1NVVZWlh544AE98sgjysnJUXZ2th599FGNGDEiePZ2d5CcAQBG+rr6dXKHsPDn7Nu3T+PHjw++LisrkyTNmjVL5eXl2rZtmyTpuuuuC5n3b//2byoqKpIkPffcc0pOTtaMGTN06tQpTZgwQevWrVNSUlK34yA5AwCMFIvbdxYVFcn+jqz+Xe+dk5aWppUrV2rlypVhH/+csNacV61apWuuuSbYt7/xxhv12muvRXxwAABwvrCS8+WXX65nn31W+/bt0759+/QXf/EXuv3223Xw4MGeig8AkKDsKIx4FVZbe+rUqSGvly5dqlWrVmnv3r0aPnx4VAMDACS2RH4qVcRrzn6/X//8z/+skydP6sYbb+x0P5/PJ5/PF3z97Tu1AACAUGFf5/z+++/re9/7nlwulx588EFt2bJFV199daf7e73ekDuz5OXlOQoYAJAgErivHXZyvvLKK1VfX6+9e/fqb//2bzVr1ix9+OGHne6/ePFitbS0BEdjY6OjgAEACcLprTsTqa2dmpqqP/uzP5MkjR49WrW1tfr1r3+tf/qnf7rg/tG6wTgAILHE4pGRpnB8+07btkPWlAEAgDNhVc4///nPVVJSory8PJ04cUKbNm3Szp079frrr/dUfACABMXZ2t30pz/9Sffee6+OHTumrKwsXXPNNXr99dc1adKknooPAJConK4bJ0pyXrNmTU/FAQAAzuLe2gAAIyXyCWEkZwCAmZxeqxzHydnx2doAACC6qJwBAEbibG0AAEwUx61pJ2hrAwBgGCpnAICRaGsDAGCaBD5bm+QMADCUdXY4mR+fWHMGAMAwVM4AADPR1gYAwDAJnJxpawMAYJjYVc5Wn6+HqfoYHNtZdrL5MZ4akBLrELpUs+KFWIfQLYVlD8Y6hC6l/nf8noADA/HISAAAzJLIT6Uyv/QCACDBUDkDAMyUwCeEkZwBAGZK4DVn2toAABiGyhkAYCTL/no4mR+vSM4AADOx5gwAgGFYcwYAAKagcgYAmCmB29pUzgAAM9lRGGHavXu3pk6dKo/HI8uytHXr1tCQbFvl5eXyeDxKT09XUVGRDh48GLKPz+fT/Pnz1b9/f/Xr10/Tpk3TZ599FlYcJGcAAM46efKkrr32WlVVVV3w/eXLl2vFihWqqqpSbW2tcnNzNWnSJJ04cSK4T2lpqbZs2aJNmzZpz549amtr05QpU+T3+7sdB21tAICZYtDWLikpUUlJyYU/zrZVWVmpJUuWaPr06ZKk9evXy+12a+PGjZo9e7ZaWlq0Zs0avfTSS5o4caIkacOGDcrLy9OOHTt06623disOKmcAgJnOna3tZEhqbW0NGT6fL6JwGhoa1NTUpOLi4uA2l8ulcePGqaamRpJUV1enjo6OkH08Ho8KCgqC+3QHyRkAcFHLy8tTVlZWcHi93og+p6mpSZLkdrtDtrvd7uB7TU1NSk1N1aWXXtrpPt1BWxsAYKRo3SGssbFRmZmZwe0ul8tZXFbo9dO2bZ+37du6s883OaqcvV6vLMtSaWmpk48BAOB8UTpbOzMzM2REmpxzc3Ml6bwKuLm5OVhN5+bmqr29XcePH+90n+6IODnX1tZq9erVuuaaayL9CAAA4kZ+fr5yc3NVXV0d3Nbe3q5du3apsLBQkjRq1CilpKSE7HPs2DF98MEHwX26I6K2dltbm+655x69+OKLeuaZZyL5CAAAjNPW1qZPPvkk+LqhoUH19fXKzs7WoEGDVFpaqoqKCg0dOlRDhw5VRUWF+vbtq5kzZ0qSsrKy9MADD+iRRx5RTk6OsrOz9eijj2rEiBHBs7e7I6LkPHfuXE2ePFkTJ07sMjn7fL6QM+NaW1sjOSQAIMFYcrjmHMGcffv2afz48cHXZWVlkqRZs2Zp3bp1WrRokU6dOqU5c+bo+PHjGjNmjLZv366MjIzgnOeee07JycmaMWOGTp06pQkTJmjdunVKSkrqdhxhJ+dNmzZp//79qq2t7db+Xq9XTz31VLiHAQAkuhg8+KKoqEi23flfBJZlqby8XOXl5Z3uk5aWppUrV2rlypVhH/+csNacGxsbtWDBAm3YsEFpaWndmrN48WK1tLQER2NjY0SBAgCQKMKqnOvq6tTc3KxRo0YFt/n9fu3evVtVVVXy+Xznle0ul8vxaesAgASUwA++CCs5T5gwQe+//37Itp/97GcaNmyYHnvssbD66QAAfCeSc/dkZGSooKAgZFu/fv2Uk5Nz3nYAABAZ7hAGADBStO4QFo8cJ+edO3dGIQwAAL4lgdvaPPgCAADD0NYGAJgpgStnkjMAwEiJvOZMWxsAAMNQOQMAzBSD23eaguQMADATa84AAJiFNWcAAGAMKmcAgJloawMAYBiHbe14Ts60tQEAMAyVMwDATLS1AQAwDMk5BuyApEDMDt+lgMGxnWWd7oh1CF1KbzY/xsKHH4x1CN2S/qX536XVfibWIXTNjoPf2AGDY4yH7+8iQOUMADAS1zkDAABjkJwBADAMbW0AgJk4IQwAALMk8pozyRkAYK44TrBOsOYMAIBhqJwBAGZizRkAALMk8pozbW0AAAxD5QwAMBNtbQAAzEJbGwAAGIPKGQBgJtraAAAYJoGTM21tAAAknTlzRk888YTy8/OVnp6uIUOG6Omnn1YgEAjuY9u2ysvL5fF4lJ6erqKiIh08eDDqsYSVnMvLy2VZVsjIzc2NelAAAJw7IczJCMeyZcv0wgsvqKqqSh999JGWL1+uX/7yl1q5cmVwn+XLl2vFihWqqqpSbW2tcnNzNWnSJJ04cSKqP3vYbe3hw4drx44dwddJSUlRDQgAAEm93tb+wx/+oNtvv12TJ0+WJF1xxRV6+eWXtW/fvq8/zrZVWVmpJUuWaPr06ZKk9evXy+12a+PGjZo9e7aDYEOF3dZOTk5Wbm5ucFx22WVRCwYAgCA7CkNSa2tryPD5fBc83NixY/Xmm2/q8OHDkqR3331Xe/bs0W233SZJamhoUFNTk4qLi4NzXC6Xxo0bp5qamqj+6GEn548//lgej0f5+fm666679Omnn0Y1IAAAoikvL09ZWVnB4fV6L7jfY489prvvvlvDhg1TSkqKRo4cqdLSUt19992SpKamJkmS2+0Omed2u4PvRUtYbe0xY8boN7/5jX7wgx/oT3/6k5555hkVFhbq4MGDysnJueAcn88X8ldKa2urs4gBAAkhWjchaWxsVGZmZnC7y+W64P6bN2/Whg0btHHjRg0fPlz19fUqLS2Vx+PRrFmz/vdzLStknm3b521zKqzkXFJSEvzfI0aM0I033qjvf//7Wr9+vcrKyi44x+v16qmnnnIWJQAg8URpzTkzMzMkOXdm4cKFevzxx3XXXXdJ+jrPHTlyRF6vV7NmzQqeAN3U1KSBAwcG5zU3N59XTTvl6FKqfv36acSIEfr444873Wfx4sVqaWkJjsbGRieHBACgR3z11Vfq0yc0LSYlJQUvpcrPz1dubq6qq6uD77e3t2vXrl0qLCyMaiyObkLi8/n00Ucf6eabb+50H5fL1WkLAQCAzvT2vbWnTp2qpUuXatCgQRo+fLgOHDigFStW6P777//68yxLpaWlqqio0NChQzV06FBVVFSob9++mjlzZuSBXkBYyfnRRx/V1KlTNWjQIDU3N+uZZ55Ra2trSC8eAICo6OVLqVauXKknn3xSc+bMUXNzszwej2bPnq1f/OIXwX0WLVqkU6dOac6cOTp+/LjGjBmj7du3KyMjw0Gg5wsrOX/22We6++679eWXX+qyyy7TD3/4Q+3du1eDBw+OalAAAPS2jIwMVVZWqrKystN9LMtSeXm5ysvLezSWsJLzpk2beioOAABCJfC9tXnwBQDASNbZ4WR+vOLBFwAAGIbKGQBgJtraAACYpbcvpTIJyRkAYKYErpxZcwYAwDBUzgAAc8Vx9esEyRkAYKREXnOmrQ0AgGGonAEAZkrgE8JIzgAAI9HWBgAAxqByBgCYibY2AABmSeS2dsySs93eIdsy+JkhKR2xjqBLVrv5MbqaTsQ6hC6l/r84+Rv1TCDWEXTJOt0e6xC6duZMrCPoku33xzqETtm2ubFdTOLktxIAIOHQ1gYAwDAkZwAAzJLIa85cSgUAgGGonAEAZqKtDQCAWSzblmVHnmGdzI012toAABiGyhkAYCba2gAAmIWztQEAgDGonAEAZqKtDQCAWWhrAwAAY1A5AwDMRFsbAACz0NYOw+eff66f/OQnysnJUd++fXXdddeprq6uJ2IDACQyOwojToVVOR8/flw33XSTxo8fr9dee00DBgzQH//4R11yySU9FB4AAIknrOS8bNky5eXlae3atcFtV1xxRbRjAgBAUny3pp0Iq629bds2jR49WnfeeacGDBigkSNH6sUXX+yp2AAAicy2nY8wdbV0a9u2ysvL5fF4lJ6erqKiIh08eDCaP7WkMJPzp59+qlWrVmno0KF644039OCDD+qhhx7Sb37zm07n+Hw+tba2hgwAAExzbuk2JSVFr732mj788EP96le/Clm6Xb58uVasWKGqqirV1tYqNzdXkyZN0okTJ6IaS1ht7UAgoNGjR6uiokKSNHLkSB08eFCrVq3Sfffdd8E5Xq9XTz31lPNIAQAJpbfP1u5q6da2bVVWVmrJkiWaPn26JGn9+vVyu93auHGjZs+eHXmw3xJW5Txw4EBdffXVIduuuuoqHT16tNM5ixcvVktLS3A0NjZGFikAILFE6Wztb3dvfT7fBQ/X1dJtQ0ODmpqaVFxcHNzmcrk0btw41dTURPVHDys533TTTTp06FDItsOHD2vw4MGdznG5XMrMzAwZAAD0lry8PGVlZQWH1+u94H5dLd02NTVJktxud8g8t9sdfC9awmprP/zwwyosLFRFRYVmzJihd955R6tXr9bq1aujGhQAAFbg6+FkviQ1NjaGFIYul+uC+3d36dayrJB5tm2ft82psCrnG264QVu2bNHLL7+sgoIC/f3f/70qKyt1zz33RDUoAACi1db+dve2s+Tc1dJtbm6uJJ1XJTc3N59XTTsV9u07p0yZoilTpkQ1CAAAYq2rpdv8/Hzl5uaqurpaI0eOlCS1t7dr165dWrZsWVRj4d7aAAAj9fbZ2l0t3VqWpdLSUlVUVGjo0KEaOnSoKioq1LdvX82cOTPyQC+A5AwAMFOENxIJmR+Gc0u3ixcv1tNPP638/Pzzlm4XLVqkU6dOac6cOTp+/LjGjBmj7du3KyMjI/I4L4DkDAAwUiyeStXV0q1lWSovL1d5eXnkgXVD2E+lAgAAPYvKGQBgJqePfYzjh2aQnAEARopFW9sUtLUBADAMlTMAwEy9fLa2SUjOAAAj0dYGAADGoHIGAJiJs7UBADALbW0AAGAMKmcAgJkC9tfDyfw4FbvkHHD4FO2e1tEe6wi6dsr8xofV3hHrELpkpabEOoTuiYfLQk77Yh1Bl2yf+f/fts+ciXUInbJtfy8eTKw5AwBgEksO15yjFknvM7/0AgAgwVA5AwDMxB3CAAAwC5dSAQAAY1A5AwDMxNnaAACYxbJtWQ7WjZ3MjTXa2gAAGIbKGQBgpsDZ4WR+nCI5AwCMRFsbAAAYg8oZAGAmztYGAMAw3CEMAACzcIcwAABgDCpnAICZEritHVblfMUVV8iyrPPG3Llzeyo+AECCsgLOR7wKq3Kura2V3+8Pvv7ggw80adIk3XnnnVEPDACARBVWcr7ssstCXj/77LP6/ve/r3HjxkU1KAAAErmtHfGac3t7uzZs2KCysjJZltXpfj6fTz6fL/i6tbU10kMCABJJAl/nHPHZ2lu3btX//M//6Kc//el37uf1epWVlRUceXl5kR4SAICEEHFyXrNmjUpKSuTxeL5zv8WLF6ulpSU4GhsbIz0kACCBnLu3tpMRryJqax85ckQ7duzQK6+80uW+LpdLLpcrksMAABJZAq85R1Q5r127VgMGDNDkyZOjHQ8AAEbwer2yLEulpaXBbbZtq7y8XB6PR+np6SoqKtLBgwejfuywk3MgENDatWs1a9YsJSdzDxMAQA+x9b/PdI5kOCica2trtXr1al1zzTUh25cvX64VK1aoqqpKtbW1ys3N1aRJk3TixInID3YBYSfnHTt26OjRo7r//vujGggAAN8UqzXntrY23XPPPXrxxRd16aWXBrfbtq3KykotWbJE06dPV0FBgdavX6+vvvpKGzdujNaPLSmC5FxcXCzbtvWDH/wgqoEAABDC1v+uO0c0vv6Y1tbWkPHNy3svZO7cuZo8ebImTpwYsr2hoUFNTU0qLi4ObnO5XBo3bpxqamqi+qPz4AsAwEUtLy8v5JJer9fb6b6bNm3S/v37L7hPU1OTJMntdodsd7vdwfeihUVjAICZonS2dmNjozIzM4ObO7uCqLGxUQsWLND27duVlpbW6cd++8Zbtm1/5824IkFyBgCYKSDJSc47++CLzMzMkOTcmbq6OjU3N2vUqFHBbX6/X7t371ZVVZUOHTok6esKeuDAgcF9mpubz6umnaKtDQCApAkTJuj9999XfX19cIwePVr33HOP6uvrNWTIEOXm5qq6ujo4p729Xbt27VJhYWFUY6FyBgAYyeldvsKdm5GRoYKCgpBt/fr1U05OTnB7aWmpKioqNHToUA0dOlQVFRXq27evZs6cGXGcF0JyBgCYycA7hC1atEinTp3SnDlzdPz4cY0ZM0bbt29XRkZGVI9DcgYAoBM7d+4MeW1ZlsrLy1VeXt6jxyU5AwDMZGDl3FtIzgAAMyVwcuZsbQAADEPlDAAwU5Suc45HJGcAgJF6+1Iqk5CcAQBmSuA155gl58BpnwKWuT2HePiLy/Kb+/2dYyXFwWkNp5NiHUH3RPnevT3CMv/f2+7iiUQmsDvOxDqETtm2ubFdTKicAQBmCtiS5aBQCphfZHWG5AwAMFMCt7XN70EBAJBgqJwBAIZyWDkrfitnkjMAwEy0tQEAgCmonAEAZgrYctSa5mxtAACizA58PZzMj1O0tQEAMAyVMwDATAl8QhjJGQBgJtacAQAwTAJXzqw5AwBgGCpnAICZbDmsnKMWSa8jOQMAzERbGwAAmCKs5HzmzBk98cQTys/PV3p6uoYMGaKnn35agUD8XugNADBUIOB8xKmw2trLli3TCy+8oPXr12v48OHat2+ffvaznykrK0sLFizoqRgBAIkogdvaYSXnP/zhD7r99ts1efJkSdIVV1yhl19+Wfv27euR4AAASERhtbXHjh2rN998U4cPH5Ykvfvuu9qzZ49uu+22Tuf4fD61traGDAAAunSucnYy4lRYlfNjjz2mlpYWDRs2TElJSfL7/Vq6dKnuvvvuTud4vV499dRTjgMFACSYBL5DWFiV8+bNm7VhwwZt3LhR+/fv1/r16/UP//APWr9+fadzFi9erJaWluBobGx0HDQAABezsCrnhQsX6vHHH9ddd90lSRoxYoSOHDkir9erWbNmXXCOy+WSy+VyHikAIKHYdkC2g8c+Opkba2El56+++kp9+oQW20lJSVxKBQCIPtt21ppOlDXnqVOnaunSpRo0aJCGDx+uAwcOaMWKFbr//vt7Kj4AQKKyHa45J0pyXrlypZ588knNmTNHzc3N8ng8mj17tn7xi1/0VHwAACScsJJzRkaGKisrVVlZ2UPhAABwViAgWQ6WTeN4zZl7awMAzNTL1zl7vV7dcMMNysjI0IABA3THHXfo0KFD3wrJVnl5uTwej9LT01VUVKSDBw9G86eWRHIGAECStGvXLs2dO1d79+5VdXW1zpw5o+LiYp08eTK4z/Lly7VixQpVVVWptrZWubm5mjRpkk6cOBHVWHhkJADASHYgINtBWzvcS6lef/31kNdr167VgAEDVFdXp1tuuUW2bauyslJLlizR9OnTJUnr16+X2+3Wxo0bNXv27Ihj/TYqZwCAmWJ8+86WlhZJUnZ2tiSpoaFBTU1NKi4uDu7jcrk0btw41dTUODrWt1E5AwAuat9+pkN3bo5l27bKyso0duxYFRQUSJKampokSW63O2Rft9utI0eORDFiKmcAgKkCtvMhKS8vT1lZWcHh9Xq7PPS8efP03nvv6eWXXz7vPcuyQl7btn3eNqeonAEAZrJtSU4upfo6OTc2NiozMzO4uauqef78+dq2bZt2796tyy+/PLg9NzdX0tcV9MCBA4Pbm5ubz6umnaJyBgBc1DIzM0NGZ8nZtm3NmzdPr7zyit566y3l5+eHvJ+fn6/c3FxVV1cHt7W3t2vXrl0qLCyMasxUzgAAI9kBW7YV+UlddpgnhM2dO1cbN27Uq6++qoyMjOAac1ZWltLT02VZlkpLS1VRUaGhQ4dq6NChqqioUN++fTVz5syI47wQkjMAwEx2QM7a2uHNXbVqlSSpqKgoZPvatWv105/+VJK0aNEinTp1SnPmzNHx48c1ZswYbd++XRkZGZHHeQEkZwCAkXq7cu7O/pZlqby8XOXl5RFG1T2sOQMAYJher5zP/WVyRh2OngTW0yzb/L9bHPxB2Wvi4XtUICnWEXRPlC/V6BGW+f/ett0e6xC6FLA7Yh1Cp86cjS3cqjSyY/kcPbzijMz9HrvS68n53P1H9+j3vX3o8PhiHUA3xEOMAC5KJ06cUFZWVo98dmpqqnJzc7WnyXmeyM3NVWpqahSi6l2W3Rt//nxDIBDQF198oYyMjKhctN3a2qq8vLzzrmNDePgeo4PvMXr4LqMj2t+jbds6ceKEPB6P+vTpuU7J6dOn1d7uvMuRmpqqtLS0KETUu3q9cu7Tp0/IRd3Rcu76NTjD9xgdfI/Rw3cZHdH8HnuqYv6mtLS0uEyq0WL+AhEAAAmG5AwAgGHiPjm7XC793d/9XZf3SsV343uMDr7H6OG7jA6+x/jU6yeEAQCA7xb3lTMAABcbkjMAAIYhOQMAYBiSMwAAhon75Pz8888rPz9faWlpGjVqlN5+++1YhxRXvF6vbrjhBmVkZGjAgAG64447dOjQoViHFfe8Xm/w2a8Iz+eff66f/OQnysnJUd++fXXdddeprq4u1mHFlTNnzuiJJ55Qfn6+0tPTNWTIED399NMKBBw8fhG9Kq6T8+bNm1VaWqolS5bowIEDuvnmm1VSUqKjR4/GOrS4sWvXLs2dO1d79+5VdXW1zpw5o+LiYp08eTLWocWt2tparV69Wtdcc02sQ4k7x48f10033aSUlBS99tpr+vDDD/WrX/1Kl1xySaxDiyvLli3TCy+8oKqqKn300Udavny5fvnLX2rlypWxDg3dFNeXUo0ZM0bXX3998AHZknTVVVfpjjvukNfrjWFk8eu//uu/NGDAAO3atUu33HJLrMOJO21tbbr++uv1/PPP65lnntF1112nysrKWIcVNx5//HH9+7//Ox0wh6ZMmSK32601a9YEt/3VX/2V+vbtq5deeimGkaG74rZybm9vV11dnYqLi0O2FxcXq6amJkZRxb+WlhZJUnZ2dowjiU9z587V5MmTNXHixFiHEpe2bdum0aNH684779SAAQM0cuRIvfjii7EOK+6MHTtWb775pg4fPixJevfdd7Vnzx7ddtttMY4M3dXrD76Ili+//FJ+v19utztku9vtVlNTU4yiim+2bausrExjx45VQUFBrMOJO5s2bdL+/ftVW1sb61Di1qeffqpVq1aprKxMP//5z/XOO+/ooYceksvl0n333Rfr8OLGY489ppaWFg0bNkxJSUny+/1aunSp7r777liHhm6K2+R8zrcfO2nbdlQeRZmI5s2bp/fee0979uyJdShxp7GxUQsWLND27dsT+kk6TgUCAY0ePVoVFRWSpJEjR+rgwYNatWoVyTkMmzdv1oYNG7Rx40YNHz5c9fX1Ki0tlcfj0axZs2IdHrohbpNz//79lZSUdF6V3NzcfF41ja7Nnz9f27Zt0+7du3vkkZ4Xu7q6OjU3N2vUqFHBbX6/X7t371ZVVZV8Pp+SkpJiGGF8GDhwoK6++uqQbVdddZV++9vfxiii+LRw4UI9/vjjuuuuuyRJI0aM0JEjR+T1eknOcSJu15xTU1M1atQoVVdXh2yvrq5WYWFhjKKKP7Zta968eXrllVf01ltvKT8/P9YhxaUJEybo/fffV319fXCMHj1a99xzj+rr60nM3XTTTTeddynf4cOHNXjw4BhFFJ+++uor9ekT+us9KSmJS6niSNxWzpJUVlame++9V6NHj9aNN96o1atX6+jRo3rwwQdjHVrcmDt3rjZu3KhXX31VGRkZwU5EVlaW0tPTYxxd/MjIyDhvnb5fv37Kyclh/T4MDz/8sAoLC1VRUaEZM2bonXfe0erVq7V69epYhxZXpk6dqqVLl2rQoEEaPny4Dhw4oBUrVuj++++PdWjoLjvO/eM//qM9ePBgOzU11b7++uvtXbt2xTqkuCLpgmPt2rWxDi3ujRs3zl6wYEGsw4g7v/vd7+yCggLb5XLZw4YNs1evXh3rkOJOa2urvWDBAnvQoEF2WlqaPWTIEHvJkiW2z+eLdWjopri+zhkAgItR3K45AwBwsSI5AwBgGJIzAACGITkDAGAYkjMAAIYhOQMAYBiSMwAAhiE5AwBgGJIzAACGITkDAGAYkjMAAIYhOQMAYJj/D8LRstY4dAveAAAAAElFTkSuQmCC", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAGdCAYAAAAotLvzAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAK2lJREFUeJzt3X90VPWd//HXEMgEMIkGzK/DEOMurkrwV+CgEQUEYlPAUlnxVy1UukeXQEkjtSLtmlrIrO4W2QNrung4/CzC7llRtioS6gJlWVaIUgFdwJVCqMSsfjEhiBOYud8/NLOO/MhM7iT3c7nPxzmfc5ibe+e+M4fknff787n3+izLsgQAAIzRzekAAABALJIzAACGITkDAGAYkjMAAIYhOQMAYBiSMwAAhiE5AwBgGJIzAACG6d7VJ4xEIvroo4+Unp4un8/X1acHANhgWZZOnDih/Px8devWefXdF198odbWVtvvk5qaqrS0tCRE1LW6PDl/9NFHCgQCXX1aAEAS1dfXq1+/fp3y3l988YUKCy5RQ2PY9nvl5ubq0KFDrkvQXZ6c09PTJUnXPvhzpaSa+2G1Zppf1Z/u5XQE7Yv4zb87rJXidATx8dn/PdXpuoXM/7np8bnTEbQvtcncn5tw6xd67ze/jP4u7wytra1qaAzrUF2BMtI7Xp03n4iosPiwWltbSc7taWtlp6SmGZ2cU/zm/5KJmPvx/R+Sc9K4ITmnuGCqKsUNn2Oq+T83XTEtmZHezVZydrMuT84AAMQjbEUUtvF3StiKJC+YLkZyBgAYKSJLEXU8O9s51mkkZwCAkSKKyE7ta+9oZ3mzmQ8AgMGonAEARgpblsJWx1vTdo51GskZAGAkL88509YGAMAwVM4AACNFZCns0cqZ5AwAMBJtbQAAYAwqZwCAkVitDQCAYSJfDTvHuxVtbQAADNOh5Pz888+rsLBQaWlpKi4u1u9///tkxwUA8LjwV6u17Qy3Sjg5r127VhUVFZozZ47eeecd3XbbbSorK9ORI0c6Iz4AgEeFLfvDrRJOzvPnz9fUqVP1wx/+UNdcc40WLFigQCCgmpqazogPAOBRkSQMt0ooObe2tqqurk6lpaUx20tLS7V9+/ZzHhMKhdTc3BwzAADA+SWUnD/55BOFw2Hl5OTEbM/JyVFDQ8M5jwkGg8rMzIyOQCDQ8WgBAJ4RkU9hGyMin9PfQod1aEGYzxf7DVuWdda2NrNnz1ZTU1N01NfXd+SUAACPiVj2h1sldJ1z3759lZKSclaV3NjYeFY13cbv98vv93c8QgAAPCahyjk1NVXFxcWqra2N2V5bW6uSkpKkBgYA8DY7Le224VYJ3yGssrJSDz30kAYPHqxbbrlFixcv1pEjR/Too492RnwAAI+ym2A9lZzvvfdeffrpp3r66ad17NgxFRUV6bXXXlNBQUFnxAcAgOd06N7a06ZN07Rp05IdCwAAURHLp4jV8erXzrFO48EXAAAjebmtzYMvAAAwDJUzAMBIYXVT2EYNGU5iLF2N5AwAMJJlc87ZYs4ZAIDkYs4ZAAAYg8oZAGCksNVNYcvGnLNX7q0NAEBXiciniI0Gb0Tuzc60tQEAMAyVMwDASF5eEEZyBgAYyf6cM21tAACQJFTOAAAjfbkgzMaDL2hrJ64106cUv7kf3Bd9zW+HnLn0jNMhtCvlktNOh9CulJSI0yHEJRw2v9EVaunhdAjtCn/mhprE3N+N4VDXxRaxeftOVmsDAICkccOfkAAAD/LygjCSMwDASBF18+xNSEjOAAAjhS2fwjaeLGXnWKcx5wwAgGGonAEARgrbXK0dpq0NAEByRaxuithYEBZx8YIw2toAABiGyhkAYCTa2gAAGCYieyuu3XHvv3OjrQ0AgGGonAEARrJ/ExL31p8kZwCAkezfvtO9ydm9kQMAcJGicgYAGMnLz3OmcgYAGKmtrW1nJCIYDGrIkCFKT09Xdna2JkyYoP3798fsM2XKFPl8vphx8803x+wTCoU0Y8YM9e3bV71799Zdd92lo0ePJhRLwsl569atGj9+vPLz8+Xz+fTyyy8n+hYAALSr7TpnOyMRW7ZsUXl5uXbs2KHa2lqdOXNGpaWlOnnyZMx+3/rWt3Ts2LHoeO2112K+XlFRoXXr1mnNmjXatm2bWlpaNG7cOIXD4bhjSbitffLkSV1//fX6wQ9+oIkTJyZ6OAAARtqwYUPM66VLlyo7O1t1dXW6/fbbo9v9fr9yc3PP+R5NTU1asmSJVq5cqdGjR0uSVq1apUAgoE2bNunOO++MK5aEk3NZWZnKysoSPQwAgIRELJ8idm5CYvORkU1NTZKkrKysmO2bN29Wdna2Lr30Ug0fPlzz5s1Tdna2JKmurk6nT59WaWlpdP/8/HwVFRVp+/btnZecExUKhRQKhaKvm5ubO/uUAICLQMTm7TvbrnP+Zt7x+/3y+/0XPNayLFVWVmrYsGEqKiqKbi8rK9M999yjgoICHTp0SD//+c91xx13qK6uTn6/Xw0NDUpNTdVll10W8345OTlqaGiIO/ZOXxAWDAaVmZkZHYFAoLNPCQBAVCAQiMlDwWCw3WOmT5+ud999Vy+++GLM9nvvvVdjx45VUVGRxo8fr9dff10HDhzQq6++esH3syxLPl/8lXynV86zZ89WZWVl9HVzczMJGgDQLvuPjPzy2Pr6emVkZES3t1c1z5gxQ+vXr9fWrVvVr1+/C+6bl5engoICHTx4UJKUm5ur1tZWHT9+PKZ6bmxsVElJSdyxd3rl7Pf7lZGRETMAAGhPWD7bQ9JZOeh8ydmyLE2fPl0vvfSS3nzzTRUWFrYb46effqr6+nrl5eVJkoqLi9WjRw/V1tZG9zl27Jj27t2bUHLmJiQAAEgqLy/X6tWr9corryg9PT06R5yZmamePXuqpaVFVVVVmjhxovLy8vTHP/5RTz75pPr27avvfve70X2nTp2qxx57TH369FFWVpZmzZqlQYMGRVdvxyPh5NzS0qIPPvgg+vrQoUPavXu3srKy1L9//0TfDgCAc0pWWzteNTU1kqQRI0bEbF+6dKmmTJmilJQU7dmzRytWrNBnn32mvLw8jRw5UmvXrlV6enp0/+eee07du3fXpEmTdOrUKY0aNUrLli1TSkpK3LEknJx37dqlkSNHRl+3zSdPnjxZy5YtS/TtAAA4p7AUbU139PhEWJZ1wa/37NlTb7zxRrvvk5aWpoULF2rhwoUJRvB/Ek7OI0aMaPcbAAAAHcecMwDASF3d1jYJyRkAYCQvP8+Z5AwAMJJl85GRFo+MBAAAyULlDAAwEm1tAAAM4/RTqZzk3j8rAAC4SFE5AwCMFLb5yEg7xzqN5AwAMBJtbQAAYAwqZwCAkSLqpoiNGtLOsU4jOQMAjBS2fArbaE3bOdZp7v2zAgCAi5RjlfPpXlIkzamzt+/MpWecDqFdl1x+0ukQ2lVw2XGnQ2hXlt/8z1GS/l+ot9MhtOvIZ5c6HUK7TugSp0No1+lTPZwO4bzC8T+S2DYvLwijrQ0AMJJl86lUFncIAwAgucLyKWzj4RV2jnWae/+sAADgIkXlDAAwUsSyN28csZIYTBcjOQMAjBSxOeds51inuTdyAAAuUlTOAAAjReRTxMaiLjvHOo3kDAAwEncIAwAAxqByBgAYycsLwkjOAAAjRWTz9p0unnN2758VAABcpKicAQBGsmyu1rZcXDmTnAEARuKpVAAAGMbLC8LcGzkAABephJJzMBjUkCFDlJ6eruzsbE2YMEH79+/vrNgAAB7W1ta2M9wqoeS8ZcsWlZeXa8eOHaqtrdWZM2dUWlqqkydPdlZ8AACPart9p53hVgnNOW/YsCHm9dKlS5Wdna26ujrdfvvtSQ0MAACvsrUgrKmpSZKUlZV13n1CoZBCoVD0dXNzs51TAgA8wsurtTu8IMyyLFVWVmrYsGEqKio6737BYFCZmZnREQgEOnpKAICHMOfcAdOnT9e7776rF1988YL7zZ49W01NTdFRX1/f0VMCAOAJHWprz5gxQ+vXr9fWrVvVr1+/C+7r9/vl9/s7FBwAwLu83NZOKDlblqUZM2Zo3bp12rx5swoLCzsrLgCAx5Gc41ReXq7Vq1frlVdeUXp6uhoaGiRJmZmZ6tmzZ6cECACA1yQ051xTU6OmpiaNGDFCeXl50bF27drOig8A4FGW7F3rbDn9DdiQcFsbAICuQFsbAADDeDk58+ALAAAMQ+UMADCSlytnkjMAwEheTs60tQEAMAyVMwDASJblk2Wj+rVzrNNIzgAAI9l9JrObn+dMWxsAAMNQOQMAjMSCMAAADNM252xnJCIYDGrIkCFKT09Xdna2JkyYoP37938jJktVVVXKz89Xz549NWLECO3bty9mn1AopBkzZqhv377q3bu37rrrLh09ejShWEjOAABI2rJli8rLy7Vjxw7V1tbqzJkzKi0t1cmTJ6P7PPvss5o/f74WLVqknTt3Kjc3V2PGjNGJEyei+1RUVGjdunVas2aNtm3bppaWFo0bN07hcDjuWGhrAwCM1NVt7Q0bNsS8Xrp0qbKzs1VXV6fbb79dlmVpwYIFmjNnju6++25J0vLly5WTk6PVq1frkUceUVNTk5YsWaKVK1dq9OjRkqRVq1YpEAho06ZNuvPOO+OKhcoZAGCkZLW1m5ubY0YoFIrr/E1NTZKkrKwsSdKhQ4fU0NCg0tLS6D5+v1/Dhw/X9u3bJUl1dXU6ffp0zD75+fkqKiqK7hMPxyrniN+S/OY+5SrlktNOh9CugsuOOx1Cu76fH/9/RqdMuqTJ6RDi8s8tmU6H0K4VKnE6hHb99+d+p0NoV9hvblMz0oVPJ7RsVs5tyTkQCMRsf+qpp1RVVdXOsZYqKys1bNgwFRUVSZIaGhokSTk5OTH75uTk6PDhw9F9UlNTddlll521T9vx8TD3fwAAAElQX1+vjIyM6Gu/v/0/0KZPn653331X27ZtO+trPl/sHwyWZZ217Zvi2efraGsDAIxkSbIsG+Or98nIyIgZ7SXnGTNmaP369fr3f/939evXL7o9NzdXks6qgBsbG6PVdG5urlpbW3X8+PHz7hMPkjMAwEhtdwizMxJhWZamT5+ul156SW+++aYKCwtjvl5YWKjc3FzV1tZGt7W2tmrLli0qKflySqe4uFg9evSI2efYsWPau3dvdJ940NYGAEBSeXm5Vq9erVdeeUXp6enRCjkzM1M9e/aUz+dTRUWFqqurNWDAAA0YMEDV1dXq1auXHnjggei+U6dO1WOPPaY+ffooKytLs2bN0qBBg6Krt+NBcgYAGKmrH3xRU1MjSRoxYkTM9qVLl2rKlCmSpMcff1ynTp3StGnTdPz4cQ0dOlQbN25Uenp6dP/nnntO3bt316RJk3Tq1CmNGjVKy5YtU0pKStyxkJwBAEaKWD75uvA6ZyuOleg+n09VVVUXXO2dlpamhQsXauHChQmd/+uYcwYAwDBUzgAAI7WturZzvFuRnAEARurqOWeT0NYGAMAwVM4AACN5uXImOQMAjNTVq7VNQnIGABjJywvCmHMGAMAwVM4AACN9WTnbmXNOYjBdjOQMADCSlxeEJdTWrqmp0XXXXRd97NYtt9yi119/vbNiAwDAkxJKzv369dPf/u3fateuXdq1a5fuuOMOfec739G+ffs6Kz4AgEdZSRhulVBbe/z48TGv582bp5qaGu3YsUMDBw5MamAAAG/zclu7w3PO4XBY//Iv/6KTJ0/qlltuOe9+oVBIoVAo+rq5ubmjpwQAwBMSvpRqz549uuSSS+T3+/Xoo49q3bp1uvbaa8+7fzAYVGZmZnQEAgFbAQMAPMLDfe2Ek/Nf/MVfaPfu3dqxY4f++q//WpMnT9Z777133v1nz56tpqam6Kivr7cVMADAI75qa3d0yEtt7dTUVP35n/+5JGnw4MHauXOn/uEf/kH/9E//dM79/X6//H6/vSgBAJ7DHcJssCwrZk4ZAADYk1Dl/OSTT6qsrEyBQEAnTpzQmjVrtHnzZm3YsKGz4gMAeBSrteP08ccf66GHHtKxY8eUmZmp6667Ths2bNCYMWM6Kz4AgFfZnTf2SnJesmRJZ8UBAAC+wr21AQBG8vKCMJIzAMBMdq9VdnFy5nnOAAAYhsoZAGAkVmsDAGAiF7em7aCtDQCAYaicAQBGoq0NAIBpPLxam+QMADCU76th53h3Ys4ZAADDUDkDAMxEWxsAAMN4ODnT1gYAwDCOVc5WypfDVCkpEadDaFeW/6TTIbRr0iVNTofQrjvzb3A6hLi88dFup0No129d8H/SDT/bYYN/N3bp720eGQkAgFm8/FQq2toAABiGyhkAYCYPLwgjOQMAzOThOWfa2gAAGIbKGQBgJJ/15bBzvFuRnAEAZmLOGQAAwzDnDAAATEHlDAAwE21tAAAM4+HkTFsbAADDUDkDAMzk4cqZ5AwAMBOrtQEAgCmonAEARvLyHcJsVc7BYFA+n08VFRVJCgcAgK9YSRgJ2rp1q8aPH6/8/Hz5fD69/PLLMV+fMmWKfD5fzLj55ptj9gmFQpoxY4b69u2r3r1766677tLRo0cTiqPDyXnnzp1avHixrrvuuo6+BQAARjl58qSuv/56LVq06Lz7fOtb39KxY8ei47XXXov5ekVFhdatW6c1a9Zo27Ztamlp0bhx4xQOh+OOo0Nt7ZaWFj344IN64YUXNHfu3I68BQAAxikrK1NZWdkF9/H7/crNzT3n15qamrRkyRKtXLlSo0ePliStWrVKgUBAmzZt0p133hlXHB2qnMvLyzV27NjoiS8kFAqpubk5ZgAA0B6f/m/euUPjq/f5Zg4KhUK24tq8ebOys7N11VVX6a/+6q/U2NgY/VpdXZ1Onz6t0tLS6Lb8/HwVFRVp+/btcZ8j4eS8Zs0avf322woGg3HtHwwGlZmZGR2BQCDRUwIAvKjtUio7Q1IgEIjJQ/Hmr3MpKyvTb37zG7355pv61a9+pZ07d+qOO+6IJvyGhgalpqbqsssuizkuJydHDQ0NcZ8nobZ2fX29Zs6cqY0bNyotLS2uY2bPnq3Kysro6+bmZhI0AKDL1NfXKyMjI/ra7/d3+L3uvffe6L+Lioo0ePBgFRQU6NVXX9Xdd9993uMsy5LPF/911wkl57q6OjU2Nqq4uDi6LRwOa+vWrVq0aJFCoZBSUlJijvH7/bY+CACARyXpDmEZGRkxyTmZ8vLyVFBQoIMHD0qScnNz1draquPHj8dUz42NjSopKYn7fRNqa48aNUp79uzR7t27o2Pw4MF68MEHtXv37rMSMwAAHebApVSJ+vTTT1VfX6+8vDxJUnFxsXr06KHa2troPseOHdPevXsTSs4JVc7p6ekqKiqK2da7d2/16dPnrO0AALhNS0uLPvjgg+jrQ4cOaffu3crKylJWVpaqqqo0ceJE5eXl6Y9//KOefPJJ9e3bV9/97nclSZmZmZo6daoee+wx9enTR1lZWZo1a5YGDRoU1yLqNtwhDABgJCfuELZr1y6NHDky+rptzdTkyZNVU1OjPXv2aMWKFfrss8+Ul5enkSNHau3atUpPT48e89xzz6l79+6aNGmSTp06pVGjRmnZsmUJdZdtJ+fNmzfbfQsAAM7mwFOpRowYIcs6/4FvvPFGu++RlpamhQsXauHChYkH8BUefAEAgGFoawMAzMTznAEAMAtPpQIAAMagcgYAmOlrt+Ds8PEuRXIGAJiJOWcAAMzCnDMAADAGlTMAwEy0tQEAMIzNtrabkzNtbQAADEPlDAAwE21tAAAMQ3Luer7wl8NU4bD5Hf//F+rtdAjt+ueWTKdDaNcbH+12OoS4uOGzdMP/STf8bJv8u9Hk2C4mVM4AACNxnTMAADAGyRkAAMPQ1gYAmIkFYQAAmMXLc84kZwCAuVycYO1gzhkAAMNQOQMAzMScMwAAZvHynDNtbQAADEPlDAAwE21tAADMQlsbAAAYg8oZAGAm2toAABjGw8mZtjYAAIZJKDlXVVXJ5/PFjNzc3M6KDQDgYW0LwuwMt0q4rT1w4EBt2rQp+jolJSWpAQEAIMnTbe2Ek3P37t2plgEAnc/DyTnhOeeDBw8qPz9fhYWFuu+++/Thhx92RlwAAHhWQpXz0KFDtWLFCl111VX6+OOPNXfuXJWUlGjfvn3q06fPOY8JhUIKhULR183NzfYiBgB4AjchiVNZWZkmTpyoQYMGafTo0Xr11VclScuXLz/vMcFgUJmZmdERCATsRQwA8AYrCcOlbF1K1bt3bw0aNEgHDx487z6zZ89WU1NTdNTX19s5JQAAFz1bNyEJhUJ6//33ddttt513H7/fL7/fb+c0AAAPoq0dp1mzZmnLli06dOiQ/uu//kt/+Zd/qebmZk2ePLmz4gMAeJWH29oJVc5Hjx7V/fffr08++USXX365br75Zu3YsUMFBQWdFR8AAJ6TUHJes2ZNZ8UBAEAsD1/nzIMvAABG8n017BzvVjz4AgAAw1A5AwDMRFsbAACzePlSKpIzAMBMHq6cmXMGAMAwVM4AAHO5uPq1g+QMADCSl+ecaWsDAGAYKmcAgJlYEAYAgFna2tp2RqK2bt2q8ePHKz8/Xz6fTy+//HLM1y3LUlVVlfLz89WzZ0+NGDFC+/bti9knFAppxowZ6tu3r3r37q277rpLR48eTSgOkjMAAF85efKkrr/+ei1atOicX3/22Wc1f/58LVq0SDt37lRubq7GjBmjEydORPepqKjQunXrtGbNGm3btk0tLS0aN26cwuFw3HHQ1gYAmMmBtnZZWZnKysrO/XaWpQULFmjOnDm6++67JUnLly9XTk6OVq9erUceeURNTU1asmSJVq5cqdGjR0uSVq1apUAgoE2bNunOO++MKw4qZwCAkZLV1m5ubo4ZoVCoQ/EcOnRIDQ0NKi0tjW7z+/0aPny4tm/fLkmqq6vT6dOnY/bJz89XUVFRdJ94OFY5dwv5lOIz95khoZYeTofQriOfXep0CO1aoRKnQ2jXb/0nnQ4hLp+19nI6hHa54f9k2AU/2/6Qub8bZXJs5xEIBGJeP/XUU6qqqkr4fRoaGiRJOTk5MdtzcnJ0+PDh6D6pqam67LLLztqn7fh40NYGAJgpSW3t+vp6ZWRkRDf7/X5bYfm+UVhalnXWtrNCiWOfr6OtDQAwk5WEISkjIyNmdDQ55+bmStJZFXBjY2O0ms7NzVVra6uOHz9+3n3iQXIGABjJiUupLqSwsFC5ubmqra2NbmttbdWWLVtUUvLlFF5xcbF69OgRs8+xY8e0d+/e6D7xoK0NAMBXWlpa9MEHH0RfHzp0SLt371ZWVpb69++viooKVVdXa8CAARowYICqq6vVq1cvPfDAA5KkzMxMTZ06VY899pj69OmjrKwszZo1S4MGDYqu3o4HyRkAYCYHLqXatWuXRo4cGX1dWVkpSZo8ebKWLVumxx9/XKdOndK0adN0/PhxDR06VBs3blR6enr0mOeee07du3fXpEmTdOrUKY0aNUrLli1TSkpK3HGQnAEARvJZlnxWx7NzR44dMWKErAsc5/P5VFVVdcHV3mlpaVq4cKEWLlyY8PnbMOcMAIBhqJwBAGby8IMvSM4AACPxPGcAAGAMKmcAgJloawMAYBba2gAAwBhUzgAAM9HWBgDALLS1E/CnP/1J3/ve99SnTx/16tVLN9xwg+rq6jojNgCAlyXpqVRulFDlfPz4cd16660aOXKkXn/9dWVnZ+t//ud/dOmll3ZSeAAAeE9CyfmZZ55RIBDQ0qVLo9uuuOKKZMcEAIAkd7em7Uiorb1+/XoNHjxY99xzj7Kzs3XjjTfqhRde6KzYAABeZln2h0sllJw//PBD1dTUaMCAAXrjjTf06KOP6kc/+pFWrFhx3mNCoZCam5tjBgAAOL+E2tqRSESDBw9WdXW1JOnGG2/Uvn37VFNTo+9///vnPCYYDOoXv/iF/UgBAJ7Cau045eXl6dprr43Zds011+jIkSPnPWb27NlqamqKjvr6+o5FCgDwFlZrx+fWW2/V/v37Y7YdOHBABQUF5z3G7/fL7/d3LDoAADwooeT84x//WCUlJaqurtakSZP01ltvafHixVq8eHFnxQcA8Chf5Mth53i3SqitPWTIEK1bt04vvviiioqK9Mtf/lILFizQgw8+2FnxAQC8irZ2/MaNG6dx48Z1RiwAAEDcWxsAYCgvr9YmOQMAzGT3RiIuvgkJyRkAYCQvV84JP5UKAAB0LipnAICZ7K64dnHlTHIGABiJtjYAADAGlTMAwEys1gYAwCy0tQEAgDGonAEAZmK1NgAAZqGtDQAAjEHlDAAwU8T6ctg53qUcS849PpdSwk6dvX3hz8z/u+WELnE6hHb99+d+p0NoV/fuBv9H/JozZ1KcDqFd4ZYeTofQru4u+NnucdLpCM6vW6gLT8acMwAAZvHJ5pxz0iLpesw5AwBgGCpnAICZuEMYAABm4VIqAABgDCpnAICZWK0NAIBZfJYln415YzvHOo22NgAAhqFyBgCYKfLVsHO8S5GcAQBGoq0NAACMQeUMADATq7UBADAMdwgDAMAs3CEMAAAYg+QMADBTW1vbzkhAVVWVfD5fzMjNzf1aOJaqqqqUn5+vnj17asSIEdq3b1+yv2tJCSbnK6644qzAfT6fysvLOyU4AIB3+SL2R6IGDhyoY8eORceePXuiX3v22Wc1f/58LVq0SDt37lRubq7GjBmjEydOJPG7/lJCc847d+5UOByOvt67d6/GjBmje+65J+mBAQDQ1bp37x5TLbexLEsLFizQnDlzdPfdd0uSli9frpycHK1evVqPPPJIUuNIqHK+/PLLlZubGx2//e1v9Wd/9mcaPnx4UoMCAKCr29qSdPDgQeXn56uwsFD33XefPvzwQ0nSoUOH1NDQoNLS0ui+fr9fw4cP1/bt25P2Lbfp8Grt1tZWrVq1SpWVlfL5fOfdLxQKKRQKRV83Nzd39JQAAC9J0nXO38w7fr9ffr//rN2HDh2qFStW6KqrrtLHH3+suXPnqqSkRPv27VNDQ4MkKScnJ+aYnJwcHT582EaQ59bhBWEvv/yyPvvsM02ZMuWC+wWDQWVmZkZHIBDo6CkBAEhYIBCIyUPBYPCc+5WVlWnixIkaNGiQRo8erVdffVXSl+3rNt8sRi3LumCB2lEdTs5LlixRWVmZ8vPzL7jf7Nmz1dTUFB319fUdPSUAwEPa7q1tZ0hSfX19TB6aPXt2XOfv3bu3Bg0apIMHD0bnodsq6DaNjY1nVdPJ0KHkfPjwYW3atEk//OEP293X7/crIyMjZgAA0K4kzTl/Mwedq6V9LqFQSO+//77y8vJUWFio3Nxc1dbWRr/e2tqqLVu2qKSkJOnfeofmnJcuXars7GyNHTs22fEAAOCIWbNmafz48erfv78aGxs1d+5cNTc3a/LkyfL5fKqoqFB1dbUGDBigAQMGqLq6Wr169dIDDzyQ9FgSTs6RSERLly7V5MmT1b07d/8EAHQSS/aeyZzgYrKjR4/q/vvv1yeffKLLL79cN998s3bs2KGCggJJ0uOPP65Tp05p2rRpOn78uIYOHaqNGzcqPT3dRpDnlnB23bRpk44cOaKHH3446cEAANCmq5/nvGbNmgu/n8+nqqoqVVVVdTimeCWcnEtLS2W5+EkfAACXsGTzqVRJi6TLcW9tAAAMw6QxAMBMPM8ZAADDRCTZub+HncVkDqOtDQCAYaicAQBG6urV2iYhOQMAzOThOWfa2gAAGIbKGQBgJg9XziRnAICZPJycaWsDAGAYKmcAgJk8fJ0zyRkAYCQupQIAwDQennN2LDmnNllKSTX5g7PTS+kap0/1cDqEdoX95v/9dybF6Qji4ws7HUH7/CHzf256nHQ6gvb5j5v7uzHcam5sFxPzf3MCALwpYkk+G38MRNz7hwTJGQBgJg+3tbmUCgAAw1A5AwAMZbNylnsrZ5IzAMBMtLUBAIApqJwBAGaKWLLVmma1NgAASWZFvhx2jncp2toAABiGyhkAYCYPLwgjOQMAzMScMwAAhvFw5cycMwAAhqFyBgCYyZLNyjlpkXQ5kjMAwEy0tQEAgCkSSs5nzpzRz372MxUWFqpnz5668sor9fTTTysSce+F3gAAQ0Ui9odLJdTWfuaZZ/TrX/9ay5cv18CBA7Vr1y794Ac/UGZmpmbOnNlZMQIAvMjDbe2EkvN//ud/6jvf+Y7Gjh0rSbriiiv04osvateuXZ0SHAAAXpRQW3vYsGH63e9+pwMHDkiS/vCHP2jbtm369re/fd5jQqGQmpubYwYAAO1qq5ztDJdKqHL+6U9/qqamJl199dVKSUlROBzWvHnzdP/995/3mGAwqF/84he2AwUAeIyH7xCWUOW8du1arVq1SqtXr9bbb7+t5cuX6+///u+1fPny8x4ze/ZsNTU1RUd9fb3toAEAuJglVDn/5Cc/0RNPPKH77rtPkjRo0CAdPnxYwWBQkydPPucxfr9ffr/ffqQAAE+xrIgsG499tHOs0xJKzp9//rm6dYsttlNSUriUCgCQfJZlrzXtlTnn8ePHa968eerfv78GDhyod955R/Pnz9fDDz/cWfEBALzKsjnn7JXkvHDhQv385z/XtGnT1NjYqPz8fD3yyCP6m7/5m86KDwAAz0koOaenp2vBggVasGBBJ4UDAMBXIhHJZ2Pa1CtzzgAAdBkPt7V58AUAAIahcgYAGMmKRGTZaGt75lIqAAC6DG1tAABgCipnAICZIpbk82blTHIGAJjJsiTZuZTKvcmZtjYAAIahcgYAGMmKWLJstLUtF1fOJGcAgJmsiOy1tbmUCgCApPJy5cycMwAAhunyyrntL5lw6xddfeqEhEM+p0NoVzjF6QjaF3HBX66WCz5HSfKFnY4gDi74uekWcjqC9oVbzf25afvd3RVV6RkrZKs1fUankxhN1+ry5HzixAlJ0nu/+WVXnxoAkCQnTpxQZmZmp7x3amqqcnNzta3hNdvvlZubq9TU1CRE1bV8Vhc35SORiD766COlp6fL57P/V3Zzc7MCgYDq6+uVkZGRhAi9ic8xOfgck4fPMjmS/TlalqUTJ04oPz9f3bp13szoF198odbWVtvvk5qaqrS0tCRE1LW6vHLu1q2b+vXrl/T3zcjI4Ac4Cfgck4PPMXn4LJMjmZ9jZ1XMX5eWlubKpJosLAgDAMAwJGcAAAzj+uTs9/v11FNPye/3Ox2Kq/E5JgefY/LwWSYHn6M7dfmCMAAAcGGur5wBALjYkJwBADAMyRkAAMOQnAEAMIzrk/Pzzz+vwsJCpaWlqbi4WL///e+dDslVgsGghgwZovT0dGVnZ2vChAnav3+/02G5XjAYlM/nU0VFhdOhuM6f/vQnfe9731OfPn3Uq1cv3XDDDaqrq3M6LFc5c+aMfvazn6mwsFA9e/bUlVdeqaefflqRiHsfoeg1rk7Oa9euVUVFhebMmaN33nlHt912m8rKynTkyBGnQ3ONLVu2qLy8XDt27FBtba3OnDmj0tJSnTx50unQXGvnzp1avHixrrvuOqdDcZ3jx4/r1ltvVY8ePfT666/rvffe069+9StdeumlTofmKs8884x+/etfa9GiRXr//ff17LPP6u/+7u+0cOFCp0NDnFx9KdXQoUN10003qaamJrrtmmuu0YQJExQMBh2MzL3+93//V9nZ2dqyZYtuv/12p8NxnZaWFt100016/vnnNXfuXN1www1asGCB02G5xhNPPKH/+I//oANm07hx45STk6MlS5ZEt02cOFG9evXSypUrHYwM8XJt5dza2qq6ujqVlpbGbC8tLdX27dsdisr9mpqaJElZWVkOR+JO5eXlGjt2rEaPHu10KK60fv16DR48WPfcc4+ys7N144036oUXXnA6LNcZNmyYfve73+nAgQOSpD/84Q/atm2bvv3tbzscGeLV5Q++SJZPPvlE4XBYOTk5MdtzcnLU0NDgUFTuZlmWKisrNWzYMBUVFTkdjuusWbNGb7/9tnbu3Ol0KK714YcfqqamRpWVlXryySf11ltv6Uc/+pH8fr++//3vOx2ea/z0pz9VU1OTrr76aqWkpCgcDmvevHm6//77nQ4NcXJtcm7zzcdOWpaVlEdRetH06dP17rvvatu2bU6H4jr19fWaOXOmNm7c6Okn6dgViUQ0ePBgVVdXS5JuvPFG7du3TzU1NSTnBKxdu1arVq3S6tWrNXDgQO3evVsVFRXKz8/X5MmTnQ4PcXBtcu7bt69SUlLOqpIbGxvPqqbRvhkzZmj9+vXaunVrpzzS82JXV1enxsZGFRcXR7eFw2Ft3bpVixYtUigUUkpKioMRukNeXp6uvfbamG3XXHON/vVf/9WhiNzpJz/5iZ544gndd999kqRBgwbp8OHDCgaDJGeXcO2cc2pqqoqLi1VbWxuzvba2ViUlJQ5F5T6WZWn69Ol66aWX9Oabb6qwsNDpkFxp1KhR2rNnj3bv3h0dgwcP1oMPPqjdu3eTmON06623nnUp34EDB1RQUOBQRO70+eefq1u32F/vKSkpXErlIq6tnCWpsrJSDz30kAYPHqxbbrlFixcv1pEjR/Too486HZprlJeXa/Xq1XrllVeUnp4e7URkZmaqZ8+eDkfnHunp6WfN0/fu3Vt9+vRh/j4BP/7xj1VSUqLq6mpNmjRJb731lhYvXqzFixc7HZqrjB8/XvPmzVP//v01cOBAvfPOO5o/f74efvhhp0NDvCyX+8d//EeroKDASk1NtW666SZry5YtTofkKpLOOZYuXep0aK43fPhwa+bMmU6H4Tr/9m//ZhUVFVl+v9+6+uqrrcWLFzsdkus0NzdbM2fOtPr372+lpaVZV155pTVnzhwrFAo5HRri5OrrnAEAuBi5ds4ZAICLFckZAADDkJwBADAMyRkAAMOQnAEAMAzJGQAAw5CcAQAwDMkZAADDkJwBADAMyRkAAMOQnAEAMAzJGQAAw/x/eGVIs3aeIToAAAAASUVORK5CYII=", "text/plain": [ "
" ] @@ -200,7 +246,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "id": "32322cc5", "metadata": {}, "outputs": [ @@ -208,46 +254,46 @@ "name": "stderr", "output_type": "stream", "text": [ - "100%|██████████| 100000/100000 [00:06<00:00, 14506.65it/s]\n" + "100%|██████████| 100000/100000 [00:07<00:00, 13426.10it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "Time: 7.548181772232056\n" + "Time: 7.5848705768585205\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ - "100%|██████████| 100000/100000 [00:00<00:00, 548627.16it/s]\n" + "100%|██████████| 100000/100000 [00:00<00:00, 440484.69it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "Time: 0.7414956092834473\n", + "Time: 0.3511643409729004\n", "65536\n", - "Module pinballrt.camera c7d3b15 load on device 'cpu' took 0.62 ms (cached)\n" + "Module pinballrt.camera 1417a13 load on device 'cpu' took 0.25 ms (cached)\n" ] } ], "source": [ - "image = model.make_image(npix=256, pixel_size=0.2*u.arcsec, lam=np.array([1., 1000.])*u.micron, incl=45.*u.degree, pa=45.*u.degree, distance=1.*u.pc, device='cpu')" + "image = model.make_image(npix=256, pixel_size=0.2*u.arcsec, channels=np.array([1., 1000.])*u.micron, incl=45.*u.degree, pa=45.*u.degree, distance=1.*u.pc, include_gas=False, device='cpu')" ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 8, "id": "837dbdd1", "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAGdCAYAAACl2fynAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWgRJREFUeJztvW2MXdV59n+tvWc8GNeeB+OX8RTHsiqiVrFFVUghVhrMm4n/AkKIAkmkCiQaJQUsWYCSkiiKqSqcIBXygZbqiSIIpCn5AiESKIkjwKmFkIhLFKBRRBQnMcVTJzxm/IIZe/a+/x/2Xmvda+21z5wZv4z3zPWTBu+z384++xzOde77vta9jIgICCGEkA6RzfYFEEIIIdOF4kUIIaRzULwIIYR0DooXIYSQzkHxIoQQ0jkoXoQQQjoHxYsQQkjnoHgRQgjpHAOzfQEzoSxLvPXWW1i8eDGMMbN9OYQQQqaJiODQoUMYHR1Flk0/juqkeL311ltYvXr1bF8GIYSQE2Tv3r0477zzpn1cJ8Vr8eLFAIAP4//DAAZn+WoIIYRMl0kcxy48677Pp0snxcumCgcwiAFD8SKEkM5Rd9WdaelnWonG7du344Mf/CAWL16MFStW4Prrr8evfvWrYJ9bbrkFxpjg75JLLgn2mZiYwJYtW7Bs2TIsWrQI1113Hd58880ZvQBCCCHzj2mJ186dO3H77bfjpZdewo4dOzA5OYlNmzbhyJEjwX4f/ehHsW/fPvf37LPPBtu3bt2Kp556Ck888QR27dqFw4cP45prrkFRFCf+igghhMx5ppU2/OEPfxg8fuSRR7BixQrs3r0bH/nIR9z6oaEhjIyMJM8xPj6Ob33rW3j88cdx5ZVXAgC+853vYPXq1fjJT36Cq6++erqvgRBCyDzjhMZ5jY+PAwCWLl0arH/hhRewYsUKvP/978dnP/tZ7N+/323bvXs3jh8/jk2bNrl1o6OjWLduHV588cUTuRxCCCHzhBkbNkQEd955Jz784Q9j3bp1bv3mzZvxyU9+EmvWrMGePXvwla98BZdffjl2796NoaEhjI2NYcGCBTjnnHOC861cuRJjY2PJ55qYmMDExIR7fPDgwZleNiGEkDnAjMXrjjvuwC9+8Qvs2rUrWH/TTTe55XXr1uGiiy7CmjVr8Mwzz+CGG25oPZ+ItLpOtm/fjnvvvXeml0oIIWSOMaO04ZYtW/CDH/wAzz///JSDy1atWoU1a9bgjTfeAACMjIzg2LFjOHDgQLDf/v37sXLlyuQ57rnnHoyPj7u/vXv3zuSyCSGEzBGmJV4igjvuuANPPvkknnvuOaxdu3bKY95++23s3bsXq1atAgBceOGFGBwcxI4dO9w++/btw2uvvYYNGzYkzzE0NIQlS5YEf4QQQuYv00ob3n777fjud7+Lp59+GosXL3Y1quHhYSxcuBCHDx/Gtm3b8IlPfAKrVq3Cb3/7W3zpS1/CsmXL8PGPf9zte+utt+Kuu+7Cueeei6VLl+Luu+/G+vXrnfuQEEII6cW0xOvhhx8GAGzcuDFY/8gjj+CWW25Bnud49dVX8dhjj+Gdd97BqlWrcNlll+F73/te0ALkwQcfxMDAAG688UYcPXoUV1xxBR599FHkeX7ir4gQQsicx4iIzPZFTJeDBw9ieHgYG/ExtocihJAOMinH8QKexvj4+IxKQZzPixBCSOegeBFCCOkcFC9CCCGdg+JFCCGkc1C8CCGEdA6KFyGEkM5B8SKEENI5KF6EEEI6B8WLEEJI56B4EUII6RwUL0IIIZ2D4kUIIaRzULwIIYR0DooXIYSQzkHxIoQQ0jkoXoQQQjoHxYsQQkjnoHgRQgjpHBQvQgghnYPiRQghpHNQvAghhHQOihchhJDOQfEihBDSOShehBBCOgfFixBCSOegeBFCCOkcFC9CCCGdg+JFCCGkc1C8CCGEdA6KFyGEkM5B8SKEENI5KF6EEEI6B8WLEEJI56B4EUII6RwUL0IIIZ2D4kUIIaRzULwIIYR0DooXIYSQzkHxIoQQ0jkoXoQQQjoHxYsQQkjnoHgRQgjpHBQvQgghnYPiRQghpHNQvAghhHQOihchhJDOQfEihBDSOShehBBCOgfFixBCSOegeBFCCOkcFC9CCCGdg+JFCCGkc1C8CCGEdA6KFyGEkM5B8SKEENI5KF6EEEI6B8WLEEJI56B4EUII6RwUL0IIIZ2D4kUIIaRzULwIIYR0DooXIYSQzjEt8dq+fTs++MEPYvHixVixYgWuv/56/OpXvwr2ERFs27YNo6OjWLhwITZu3IjXX3892GdiYgJbtmzBsmXLsGjRIlx33XV48803T/zVEEIImRdMS7x27tyJ22+/HS+99BJ27NiByclJbNq0CUeOHHH73H///XjggQfw0EMP4eWXX8bIyAiuuuoqHDp0yO2zdetWPPXUU3jiiSewa9cuHD58GNdccw2Kojh5r4wQQsicxYiIzPTgP/zhD1ixYgV27tyJj3zkIxARjI6OYuvWrfjiF78IoIqyVq5cia9//ev43Oc+h/HxcSxfvhyPP/44brrpJgDAW2+9hdWrV+PZZ5/F1VdfPeXzHjx4EMPDw9iIj2HADM708gkhhMwSk3IcL+BpjI+PY8mSJdM+/oRqXuPj4wCApUuXAgD27NmDsbExbNq0ye0zNDSESy+9FC+++CIAYPfu3Th+/Hiwz+joKNatW+f2iZmYmMDBgweDP0IIIfOXGYuXiODOO+/Ehz/8Yaxbtw4AMDY2BgBYuXJlsO/KlSvdtrGxMSxYsADnnHNO6z4x27dvx/DwsPtbvXr1TC+bEELIHGDG4nXHHXfgF7/4Bf7jP/6jsc0YEzwWkca6mF773HPPPRgfH3d/e/funellE0IImQPMSLy2bNmCH/zgB3j++edx3nnnufUjIyMA0Iig9u/f76KxkZERHDt2DAcOHGjdJ2ZoaAhLliwJ/gghhMxfpiVeIoI77rgDTz75JJ577jmsXbs22L527VqMjIxgx44dbt2xY8ewc+dObNiwAQBw4YUXYnBwMNhn3759eO2119w+hBBCSC8GprPz7bffju9+97t4+umnsXjxYhdhDQ8PY+HChTDGYOvWrbjvvvtw/vnn4/zzz8d9992Hs88+G5/5zGfcvrfeeivuuusunHvuuVi6dCnuvvturF+/HldeeeXJf4WEEELmHNMSr4cffhgAsHHjxmD9I488gltuuQUA8IUvfAFHjx7FbbfdhgMHDuDiiy/Gj3/8YyxevNjt/+CDD2JgYAA33ngjjh49iiuuuAKPPvoo8jw/sVdDCCFkXnBC47xmC47zIoSQbjOr47wIIYSQ2YDiRQghpHNQvAghhHQOihchhJDOQfEihBDSOaZllSeEnAFM0WrN0T0jMSF9Q/Eicx9jADMPkwx9ahwAQMpp7EtRJLMPxYvMK8yJDoSfzpf8yeZUCLB7PXn9sB9h6nEPKGzkNDEPf44SQgjpOoy8yNzF1oZMBpOpHFo2nXxaTBS59YpUTuh5ejPVFENA1Uh7avLgNZg4MFWRpovKTNYegRrD6IucFiheZG6i6lwmUzWvzPgv/mwaiYey5cs6P3UCFdB2rbGIKeFovbL4teQmLXSlBKlKk2khi65Hi5m9JooYOYUwbUgIIaRzMPIicxaXKjSZS+EZY3wU06/lHABio8d0o4rpPNdUTCditOhoK88b128S+4l+Gh2FSRmkYSXaRsjpgOJF5haqzhWuVqlCtXzyakenhtbr66eeputxSnyD19NIIXphM61C1qdQsf5FTiFMGxJCCOkcjLzI3KGHSSOgTrtVKcSpIxjTa7Rvm9vwZDkN28Z29Tq/vaYW52DweqLITkR8NKajMECtrw0epb3XZehE1M+nz88ojJxEKF5k7jJVncsK3HSFJkjHpZ+zb2ZSC+s1WFnKdgekJNabqPYVCU6QYlQuQlOWPpVYZt5iLz0ciSgpYOSkwbQhIYSQzsHIi8wN6pRh4DBMoU0a8T79RkG9xnb108JphinF/swlLfuUkh74ZUp7YL0iD6Ixkxi/JSKAMTBFUT1OmTkQjgsDbCQWPx8hM4PiRbpN3EVDD6pNOAwbda7MhA7FE6hVTSkuJ2KX78ceX5bt12DSYuHETqdCjRIdLYb1OUxRQNT1GFUnaxMyl0507w9TiOTEoHiRuUmLSQNAaOTQHed19w1Nv6LTbxeMmZ5/KlLnmUIgnB1eiVsQvVlRk1IVGfJQwMrSm2B62OsDYwchJwhrXoQQQjoHIy/SXVLW+BaHoU8hRtFJfMxUfQ97RUnTGVB8MjtuWOIoq5T+oz7dE1EPXK4jMhED1DWu6idv7uthRvVGVAOTTVlGbsXM1cFY/yInCsWLdI9Ut/hUvSolRHo/XQezj7Wlvl8xaq0zTXO9QqLnMP2k22IRaKlzhfuaUOREwlZYoupcdn0pAApf08pK7/vIMj8eLMucEEqG5pg41r/ICcC0ISGEkM7ByIt0C22wSG6uI4g8T3fSSLkLU30P9TFxpNTjcRwxtR7T7zZ33nhFS6QSN9xti9j0fln02C5HvRFDa3yubPYGQFE/H1qNHFX0lejKASTTl4T0guJFOktyni77xanHc2XGp8J0nSzPwlRhFqUUTVq8AoEyfSzHzETgggto+YKPhUqkKXrx8Smxc/UrtVyU7h6asoSgVFb4Ut1fUSLXTCEGXTl0/Yvd6Mk0YdqQEEJI52DkRbrBFCYNF3XZ6ECn/fI8PD5XpowsU5FXmEKUVNowjsbaorDUtSuk1QySXh3QEqSYrBl5tT4u09GXqGjLlLpJbxYcYwDvRCzURdmUIpAcCxYaOHQKkeYNMj0oXuTMp0e3+LiLRt+pQnuMThWqmpcEtTETilUkZoEQtYnPCTgPG7Q8h7OlayehJXC/iz9HfIxKGwpKGGjxiS6gFiKTZ7WA1U/kzm2aLsS4BgY0U4i00JM+YNqQEEJI52DkRc5c+kkV6n0Dw0WG1lShjsiycGyXJKIwRFGYGON/9rWYNFrTgkD4k/FkDFaOI604yrLP6SIqEx5Tqv2MVOYMoEqdws4BlqljyuqE9tqLwt3fIIWY52i4EFvHgDGFSKYHxYt0GytE1l1oU4p5lk4V6vqXFS6bKtT7Bfb6zItRhmDQcyBk9rjUckTQ/P1EJq4sBTo9aATB9YSpxBZxywCxolKqqSqjWphRogaU/jx57rpvNFKILsVbNm30hJwAFC9yZtJvnUuP5crzMPLShgvdOSOKrkRHZU7klKgZ46dBMQaSN+te1X7q+pUg9YzC4k3TjcQik0bQjSkQLPENN+IorBCYzEc7AtVk196bonT3yRR17cqesPC1rKquZZ8z85N1Fgg7cQD9178YfZEErHkRQgjpHIy8yJlHHXX1tMTX+wXuwsDqno7QGjWuIKUYpg1FR1vueNRpRLjzBZGVXYwHNrcFVCdimQeAMjo+GICsXZCSiMpQRWTqMAMfUQkyv1+9b7UePoUI1O+Nv2D7rMFA5kyq5czZD6dR/6L7kDSheJHu0dZwV6caU8YMe4yucWWZT41p27tKDQZpwlrU9BiwQJh0bUzRSB22aFbPFGOKHKHAwNQ5OUTjukxa2IwAKGHqx5Xg1D8UULr6lcmUiCALBSza5q4k6PBRmzesMBnp3UYKUHOA0cBBmjBtSAghpHMw8iJnDrE1vs2k4XbXEZVpmjR0OjBOFdpzaWNG7lOKDZNGvVilGRFGXrmKlux+qVQi7LbodfdwGwYBUts+ypSh99UzIptMwicuxF2aKBt8I21YR1Qi8AOWa6u8M9trwwZK/zxxM19RzszUAGYgmUJ0Bg7nYiSE4kXOFFLuwpY6V9jOSa23+wLhOC0lZEGNq66ROfHJsnBZ17XiNGG9n6iaV9AEN3AbqtcY0V4Li/Zr2U1vsDVCPQGkEzJtj49d6gZhy6nSr3bOQ5FwzJcd64WoBmZM4FB07wsAMSq9mBoDBvTsRG9TnPWLTN8PMm9g2pAQQkjnYORFuoGKqIxKLzbGdaUGGZswcoJKGwapwmhZVISmo6fKwFE9lMxHbmL8dYYmjuhf9DBmTNOvUZ3M/2MkTA86F2GmDRtq5mRTRz7BNdjjS+8cLODGgvmeh34+LzcGDPApQG2aKYrq/dDmi2gAM9B7GhUjhuYN4qB4kdmlRwuo/upcfjlIL6o6V8MOH4lVKm0oUarQ5ihsmtClCNV+YuD3S9S8grRiQqSm7TS0p9KdMFQtqVrhNnhbvTGNub+M2t2lGpVAGSBwHgZNewHXUkoy45yLyLOg1VTgPkwMYAZaUoiFe+JoDjAWweYzTBsSQgjpHIy8yOwxVQuofkwa8biuhEkjbg8lJoq0dKow85GXM2UEaUebUoTfL7VsowqoiCpIG6buR+/b1UbgKhT13OKjMhEDYwOVqB9ifAmSMGw0UojKgWn0vc5iF6K9IabpPrQpQO/dSKYQbQTeGMDs3gSmD+cjFC9y+umRKqw2m9AtOFWdy55T17PiOpcTpSwQpdZUYSRk3ipvAus8Mi9OkiMUOVXrSg5kRg+3YX1cK0G3DLU6chraL/hKEHT3+VrUjKlunc3MGQOjTujTiVEKUdvokflLVQ7D4EeD3ZbriSrr9TqFKAaBCxFID2Au4H/EsP41L6F4kdNPKtpq7OO/+IJuGZkBMr/sGsraqCsVeelxXpmqqeSmPdrStaxEpBVGYn6b3y/cB5g62upbyCLhErXeRGYNY2tbBureCFBYURM/TsydvHkhqSjM2eh1M9/4TLGFvvBWeeOMHRJ2os9VMUy3kQouiOO/5juseRFCCOkcjLzI6aW1EW2zkW6j6S5QRV02vagb7tpf9tphGHTVUOeOBxirfoZBqlBZ4P0klXUNzAYHapuOtqpluOObXTbUisQt6RmFQTkCo2a7EqcT7fO0pQ2zuhbmUpzGOxH1+hK+LmWyqqGv7cQxiXQKMXpG6ClWtNuxLP1klmJCd2GuI0sJu2/o6+Hg5XkHxYucHvqtcwFVSklb4qMJJF2qUItXLVap8VwSNOM1ajmLrO6qW4aqU0lmQqt8i0hBWegD8dLjxICEkPW3Laxtpdbp+haCL3IxcD8GjBG/7ETNV7d6mTnsHuH6DDJpl3ukEAE1hsv9p3o/7GSWQN2JA87I4USy0dxYmTc4/mvewbQhIYSQzsHIi5x6prLEA4G70EVdekoT61ILbPNZ012oI7EgHZhIIZooPZipyElPe6Ijrdxa5+v92lKFWbhPw3Ho7k16ObWvS8HFEZiOwnTgVCo7ezAo2acNq4fGn7Twy2JM8tdtdHT93PU9iF2INgqqu3C4lyWqT2GJZicOoJ7/S3UMkcxFVT1nX6aDY15A8SKnl+k03HXik3uRy/OkQDXm5tKdM1SH+LCRbp0O1N+NebTNHqMdhVqwDEIhc8eH37lQj6trat6aeLtOkwVdNNx/lGDV/5gSumuT6rahBMpApXERjfvyy6buo4F692wyvPxAwOx7CuVClELlVKvnSXaiNxJ24nBXUneod2P5JBij52pjMcYwdTgPYNqQEEJI52DkRU4tdcrQmTQ08UBklwKsI6ipZkU2ibm5WprsxiYNd7wemBw09o1ShZGjsNWYkVxvz1c/bhgP0stBJBbnEVNuQz2mCwinOTHwvQ1VzGSNHNYUYYzRM3XBWFciDMr628IUQcBa7WHfB9WJA5KF06iYOHqsnzPuxGHH94lU5o3Mpx7dTSyK9s4bJd2H8wGKFzk1xO7CVJ0LCAciq30qm7RKKSYb6yrxq/dptH6KlwNRqpvspmpeiVShXtbCVOpUof0Oj9KGYRsp9JU2DLa5/9SbJd5o74GE83EpIdPtnNxtK71Q2ZPpVKFNO2aT4WBmP5tXfWSimS9y45r3CsrKku9+xGRh9/m4E4dbr4qFJvPOw9z330rVv+g+nPswbUgIIaRzMPIip5SpZkVODkS2v7gTTsTWhrt51uxbqMZ2JQcfp8Zs2W0qVVg20oZxShH1c8L9HCxzqEjLBNGW3s89djcssc5uahvnpYOnoGGu+GNMvJ9djq0XqcdAOVBFXwBcBBYENJnaW4WcOm0YTqPiQ0HJ1DQuUZNflzoEQudhIeF+dvBy+5zTZI5B8SInH50CdKsSQpTp+pMWqKxpiW+bWDKacyvoW6gHH6uGu7G7sC1VWNrlvGmNb6tt6XNBr1eCpUs3gZAgLVpum1p2HTZKBAIl6rsfmfEd4o0ogYjPnBaseJ2ueQFwj7NJ+Ma60YtwFvq6ttbeib70++v6l4kfK+ehfa+16zAzTes8619zEooXOXlM0UXD0aPO5Zf1l5sXslSdq1pf7xMIkTpfW8Pd4HEUbalaVlzXCmpbkbDp/dw+LYKlLfXhvUysU/gZkhFa4/W5VXASjOVSRg5xGpsWrDZRqxoC+/sW7OYiL/iI287wrM0XbjwX2iNr+B8+sXnDTZ2ifui4+lfDwFE/EQVszsCaFyGEkM4xbfH66U9/imuvvRajo6MwxuD73/9+sP2WW26p3Uz+75JLLgn2mZiYwJYtW7Bs2TIsWrQI1113Hd58880TeiHkDMBkobuwjrpcCrBOGdrPhatz2UirXjZq3+Sf3V7XuWzvQh99ZcG57cBkmz4UPRVKVkdcdaRV5j7qstFTqZb9vtVftb+p9slNsA0ZXGQWPHbHJc45UP/Z8w/4v2A/FQH6aw6Pia/Vvu5SvU5k6jpsGjV4DnXPsnCdf33hNuT1Xx392gHk+v1pvFfWLaqjbt0xxX3G7Hb9HmfBZyP4bLVNuUM6z7Tf2SNHjuCCCy7AQw891LrPRz/6Uezbt8/9Pfvss8H2rVu34qmnnsITTzyBXbt24fDhw7jmmmtQtI2YJ2c+8dglK0Sm+eUSfLGoLyMvapEA6fFYwReT/6LUAubSgXUNSKcHbYd4J1yZ2hY81qKUEgItCF6Q9Bd/ICpZJFgpMYwEJxagVlFqO19SaBG8bi1g/tqsYEXPkXsRKwei+1Zvg7q/bnhA/R7YoQzWTKOXe/5QqZeNHftnhcv+WIp+HDU+mplRP6hMch/SPaZd89q8eTM2b97cc5+hoSGMjIwkt42Pj+Nb3/oWHn/8cVx55ZUAgO985ztYvXo1fvKTn+Dqq6+e7iURQgiZZ5ySmPqFF17AihUr8P73vx+f/exnsX//frdt9+7dOH78ODZt2uTWjY6OYt26dXjxxReT55uYmMDBgweDP3KGoCMn/Qs33kfP0ZX49RyYNHTKKPNppziV6DpnOKdhmB70UZlO+0XpryByMWEUk4pobBqxJU2YjIB6nU9tL3NBOeD/JG/5y6Q9wmp5Dhth2ShLp/zibX3dj+CeqCgsb1mOUo32PfPLPaLsPPfT4rhUYR6khuOoPUhNp9LZjL46z0l3G27evBmf/OQnsWbNGuzZswdf+cpXcPnll2P37t0YGhrC2NgYFixYgHPOOSc4buXKlRgbG0uec/v27bj33ntP9qWSE0W5BXt2i2+bWFLZoBsTS8auM7tsx3PV64NJJqNze0s9lHOw6S4MG+tWyzZd6Pfz65sOQ7/csMfb51f7Ba2iMgmch22toqod1PraNVc5CtXrzOAt8aqrRgZvMoyXpZ6FK7UtcCHa588MAHHzhyEXZWRUnegFMPb9EUDyTFn8S+Ui1K8h8b6rz5JtUCwZ/IszGVBb4934r9h9CLD7/BzjpIvXTTfd5JbXrVuHiy66CGvWrMEzzzyDG264ofU4ET9BXsw999yDO++80z0+ePAgVq9effIumhBCSKc45eO8Vq1ahTVr1uCNN94AAIyMjODYsWM4cOBAEH3t378fGzZsSJ5jaGgIQ0NDp/pSyUmiMdUJoCIilbIxGdKDl00UOcXRWrUYzJBsjQL2h74epBzPzaUGIjfGZukoxp5LRUqBe9Aek4jC7OPq+dE72rLroZaj33HB1CcGfmCywI19khIwpfHBkoEfUJz726uHVfllUz8WFzmFQ6JMtMJAxIeCLlqb9BFZMIWJERh1T01uoqlT1HuvO2/Eg5btvG5F0Zhx2T1nlvm+ib2w18axX53klNS8NG+//Tb27t2LVatWAQAuvPBCDA4OYseOHW6fffv24bXXXmsVL3KG0Vbn0mk7i04ZAr3rXH38hRZrqPW1kOk6V2SJb7gLa8FzjjiV3nOiEtR/jBOutjqQd9chsL2XAwKxf3qbPl7XuwYEMuiPaRyf+huIzjEDR6KtgZWp4wJXZvw4uq/G/2hw9UejrPO65pVl/TkP7XFBXTT67GXNxs+sf81Nph15HT58GL/+9a/d4z179uDnP/85li5diqVLl2Lbtm34xCc+gVWrVuG3v/0tvvSlL2HZsmX4+Mc/DgAYHh7Grbfeirvuugvnnnsuli5dirvvvhvr16937kNCCCGkF9MWr5/97Ge47LLL3GNbi7r55pvx8MMP49VXX8Vjjz2Gd955B6tWrcJll12G733ve1i8eLE75sEHH8TAwABuvPFGHD16FFdccQUeffRR5GqaA3KG0uNXajIFCISDRdX6RnoxmAVZGThik4Yq6MdTnQQmDZWaSzbS1RETVBSGMFXYMHWY5jnscmDYyCXcR6cE3bVJmJK0tyeLUlkCPx9XJr49VAmfpivhIj+gSgHa98QUfr29P+3U90qUScNeg9suKtvmU4qSe8OIiJ8nrIqqxGfqchP2PczVjVdpQsn8tCooy/7MGyj8Z6ef9CFQnZfpw05hRLr3jh08eBDDw8PYiI9hwAzO9uXML/p1GNrahKnrFK7m5JeN2q8xT5edXLI+XzA310C9X+b7F0qdLizt4wHj+xHmBlL/TAsa7qo0INB0GOrGvLYJrR+wa8+BUHzcegn3MVDC1kOwnMCF/1saMV48SriDjMDP31ULnK1zmdILiSla1peVsLltBZC5ZQmPiY7PClHbardiAWTHxR1vO9GbSUFWlDCT1pYo1SSUAExRuuNRhssoSidARi2jFL9cFBA7B1hZVNvqhgcifhll6et0pcB2LRa73L2vwk4zKcfxAp7G+Pg4lixZMu3j2ZiX9Ic2WWjzhRWilEkDCMbv2MfJCK2HRVo3XtWzIqes8UEUpQwT8UzIfrn5OLWfNkHEVnn9PGUdbQVCZiO1TNS57f2JBKulCi2BVR1BV3k/y4ipN/joU01rHK7vgbfKGwShVxwMBqaR+l6VXrSryMmerDZy2I+IiIqwMogVKONNNyb6HGgLfDS7Zmjq0JFY29QpGTj7csdp+V+FEEIIOXNh5EWmJk4V9toPCB2GrtZVH9/4JR1GTsF6NzA5qnO5H+xNa7wbjGyQtr3rupZzyenHcMeH+0XLLtWnt0kU7fn1YYTWEm0Zb6FvHaAMAKXxky5mqOzx9bkM1DYTn0iPco5WZ37ZRXUCdeP8QGSpL9W9X3pbhtAqb8+bm+p8KirzlyOu5iXIKuu83Udb5/PMpRpbBy/bOWJS8371qH+ZzHDqlI5B8SLTp63OZTdrUdJtfqqNaHTSSC3rOlcwtsvXuZw13q43UF9ouv7UIlDK9m2Pacx+DITpN6PSgPZ81piRI6xzxYYN9zgSrJR49ULPKKyyZ+HsyPZxSrAEbU9kJEoHKlOG7sIhelumxnyJqWpIqO+FmturMnPU++VKvLSQZeJubjjzMoCidIKZnLQSgMmzqnFGat4voLXzhpu1k3QGpg0JIYR0DkZepDdqMHL1uIdJQ8+9pK3xU3TSAGw6z6YJMz9gFaFhA0Haz4SpwVylEY3aT0Vk8YzGcaow7qRRPU9zOdUPMUgV5hI4D5H7SKwyadgIRGDscq9UYRQ4ibXNl0bNjFylDI1+3cEMyS0pRN29Q0dUaln7QGxK1flyVFrW6Ki09FZ5H7H6SMzlGmsbfbUYGjam7Htol+3noyiCVGPQ91AynzpMmDeCvoc0b5zxULxImul2HWhrvmsf6/P19WWk9otEKWj1ZF1uRh0H1F+WSthMc73uiGEfBw7DtlpYVL9qTyFq56G4tCFygclrd50Sr4a0BH2g9EafppPCeCGr85625hXUv7RzULv+ytqNp4XJbQNSPwBMnd50Qq8MfVps9JAAU/ehEidS/nyBO1AQpPyCFGmqdVS03o3Xcqlt8a5CXf/qd+5Ajv86Y2HakBBCSOdg5EV6o/vAAVFEpVKAeluq+W6QXvS/kgOHoR1gXEdhoWGjLdVYrbYpQxdV1c1w7X5JF2GcNmwxdgQpRNfvrzlmK0gV1nNv2eOhIi8zUCIbKOtbWCKLu2nUuClHGuvhnHFiqnG5QGVwgDFA4aNRY7tlQEUwuUom5gjSg4HbMPPr4yjMxMGOisr0en9MFUXa9LOUcBGnCNTnxXcFsVFcY9yXPXlquXYe+qhQRZyZqcZ92WNi80Zj3Fd9oeSMhOJFmihrfHpzJGB2nRa4YLmPL53gCyh8rGsltuFrtRwJXPRlmba9NwUqnndLb4uXXZ3MWr4b9S+bNlQuwlyAgRKm3pYNlMjrtGGee/EyRhqClcpYiRiUtWuvLP29LpHF44hV0tB3DzEiPtWYAcjEp+C0ezGDd+qpZTdvmMruhWlDv+zcgTbVqAXPbtMDrRvOQxOlEf3JezkPnTibzB9fZq4G17PzvGH9qwswbUgIIaRzMPIinoRJI5juBGimAC1Rq6bk/na/Xg7D+njdBspNoQHUU2ugPl7V+etf8i4iUr/yey3bhr72fK0mDZeeDI0IrdGWGteFOk2Y2cgrLzEwUOX68rz0452NeHdejY3ESjE+sycGWX3NRaEig4FE9OUMDnogsfHXJgZGD3PSqUI/5KpqABxER97V2IioVJrPGjTcmC/7tDZdiSp96KcGM4idh33PuGyXdVSlzRzKGGJKH+BZ56Ef9xVFWC4LwcHLZxIUL9Iktsa3keqkYY+Pv1j0wGK9Ta93glCJVUpwQrdhnGr0tRdRbsW4ztXonDFVnStYlkb3DATHKDu8rXHlgiwX5EqwbNpwICuRq7RhFolXacVHDIpSLQfpRfXlWgsYUH85qwkjvXD4+pfJ6l3sNYjxWTKdAlQaaWtZ+v0x6v3xneQT61PpxcwPJE46D534TMN5aNPZIhAtPi6drdS4LJPW+er+lG7QNTmzoHiRClWzarSAammmm+yk0VjOGoLlxSdTAhMLUbqTRtsMx1X9C8GXY1Cn0rb5aH2qw3ujBZQay9U0cNg6SjR+y5oyBkpkWrDyEnnml61gDWQlBrKwDmPFq1SCNVnkQOG/XD21wuT2GB8VQn1XG2WVl0yqTvTan5CKttT3vgiqY6IxXEAcecU1LgTC5Bv4IjRv2M9HKQ3DxrTMG/axO1+P+lebdZ71rzMW1rwIIYR0DkZe851U2k8/jqdVj2tYpmU5Pr9+nKitBYOS7S9sHUXpOpmuranUld03ft5G2jBaryM5xNujY4A6+6aiDj0wWl+P/WWfGQlchZkRlyrMjOCsgUkAwMiig1iQVct5fayLtsoc7xXV/67/++5ilHaCsjKD2MhAUF2Mnacqy3zKK/NpwzhSCh/7Ecf6fgR1MdOMqpLpQKPrVdJ429NWe3eH/TVnPnrUGdLgXLHz0FriddPehhdTkfpMkjMaitd8JvofNtUCyhE33gXCfaYyaaSIOm9ow0bcVT4UMr1sRU0dD7+9Ol+4jHg5sS0QLAMnVmIS+yXP7es7JquMGMYZEfxynpVYetYRAMCyBUcwOvRO8BqsnaOQDG++dw4A4Fg5gP3Fn7g9yqLKv1XiWKJ0k33qupKka1nGtIpPLOSBQUPb6LXIRfctFiidakRK2DK4GaOrcWpKcIzx6TqdDhS13m5zL1Ct1+YNiD9eCR2AyrzRS+jIGQHThoQQQjoHI6/5jjZpJLpouChKGzD0sSlXYpxmVP9Kz2gLwfrAYZgyg/SIosJUYVvkZoIIS29z56z/TUZu7jjlLbeuvUz1LKwji1xFXtakMTz0Hs4dehcAsHzBIawcHAcAnDtwOHiO/ZNLcLx2OEyUOSbqFOI7Rxc68wfKLIhujLqeapblxGswiahKv2VRtAXAddvwg4wlur/1ehW9uK4c+hyWzAQpSWeVd8epSKwf23wwS3OPKKywh5hmnKVMI0l0FEhmBYrXfETVudos8Q13oe6kEacL42Oi9a1usJZ9fGquKVKBwATCY1q39Z02RHhMfA79rz+fNM+HSrwyJyJxqlCwcPA4AOD/LDiKpQvqtOHgIYwMvgMAWD3wDs4yBd4t/f+iE+UgAOBoMYj3imp5YnIAR49XyyKC0oirrZVZicw6D0sJHJcNAbYvQYmHFv3ADl+b7sLaVvMeBOujx7ZdVHxufxH1dehanUp9BvulbPOBUCVs81Koz3xR7dNPs15ty0yJIjltMG1ICCGkczDymm8Y04iyHD26aAQmjV49DIFkRw1JmTZUqjAwW7jZl/22wEWofr2LfikqWus3BdiI1hIRXpCCjB43HIYubQgkpzNBZdI4p04VLhs6jGWDhwAAIwPjWJ5Xy8uzSazIz8a+otrvPRnHewNVhPXuggU4Wi4AABwrBvC/ZWXeKOqxVyaK+Kp7qlKIRRhRoR7rFb/WIG2I8D6FXTn87Mm9um0APvozQcQaTZXiHIXWderTiP4zJi7tF2Aj+DiNmNovThbO1HHIFOKsQPGax8R1LpMSH/s45UTslTJMpQbtY9dkVwue39enACNBQ2/hiIUoaOabEiKTPkd8bskkWh+lCt01qpqXctBpQQGAJQsmcM6CowCAcwePYPlAJVjn5oexPJsAACzLFyI3GZbnQwCACTmK92QcAPCeDOLdWrwmigEcnaxE7Y+TeZBtzTJxoqK7vUMv1+6+ZNpPOQqnqpOlamN6MspGCjFyK+pu89qxGacTpxyw3CpUCD57KBDWxYoozZhIIZpApLV6k9mAaUNCCCGdg5HXfKJOGfYcz6X3BcL+hUDk9ouObxvb1WbYSA1MBnz6LxURIT22S6JjgghLt3qKIzpFPBi5NcJTPQwlE4TGA/8L3g1SzgQGwNkLKpPG8FBo0lg+cBAAsDw/gqV55SgcNNW/Q6aKqpblk3hPqmOOyAIcKauI7N1iAY4UVRR2dHIQh97LUGq3ofFpTNeCKTPe+GCAOOL0EVXa5GGiYxouwsT7ZlOIJvE+GKTnAHPvp4uQgMYsy27Zv/GNXoeFfk/s5yWrTBtqvY7Z3H5V76r63FGkFZs3mDo8rVC85hNVUzo30Z7JSteAtCFi9n/EsoRkGdwMSaWoeF0dnxvfyTvPvcXY9iS059P/g5fwX/zi55QSqf+j+uA5TRD4c2vHmhhAOeqM+KcKO6LDdVyAmPq54I+B3w9tx5dwNSIE2/x8XAK/XJYGWQa8e6wSovGJhfjjQFWnOjs/hkV1qvAscxxnmXG3PIAck3Vh549FgT8USwAAf5hcgv3Hq+X/d2wRDrx3NgDgyLHBqvu8u49ePaSEn8Or9K8NZdUUV7v9/es2wfrgPdDf0/pxtGzi5eCxEpXk8epNtI/LtmX1OSrVcZHVXdx73zvlJ8Hztuyr11O4TjsUr3mMlOIbu5aZ+vVchlFUqWzBgeDBRzfiZ8CFSFiDqKbKVU9cRwOipsJQQmascLkvGrjivBj7H1RfUlCRlwRaWG8HYMKu5RJ8gYaCFYpkvZgZJXhWFNS1KSFwomaFEdX+epLJ8YmzMFS3hFqYH8fZ2TEAVryq6GzQvItV+UL8b1EJ21hxNsYmhwEA/3t8GH84thgA8PbEIhw8NuSfB+rS7LXCXrN9bQam9CKvX6spVZ1rClFqCFO0XL2PdlkdXz9PsE2JjAmEA6EwSWJ5KvQN0WKmlgOxapmkMuguz3rXrMOaFyGEkM7ByGu+IfpndvTbpf5lWTnB7GR+1SBl+8u0mrLd/urM1S/QvGfqxEZBYtSvXyMuUgkmJCzrn++2JGFM8hd8mPKTOtqx16nqX3FEFh0fP64eIHhOPVmjTkmGkZdPzUkJNyBWsiKIvooywzvvLQQAnJUfx8I68hrKjmPQVC960BQ4VB7DEamiqrHJ/+Mir/3HFuOPE4sAAP/vvbOrKVLgIzydrnQRn44K7XXX/5rS+LhYAJ2uDaKl+D2YKsIKIjKpo7VExCoSRmulP5dR21rR6UW7nEpTx8dYyin2nQqmDGcFitd8ROXqgvqXEjMvAnUK0U3Ulyn7tKgxOeLrW8FzmDCNmLyOetkJWf3YfokW4rOOysEsxqfzbFpMp69EfSFK4gvZCZ7WcvXdmvoSbgqeCbdZkS50t/cqhVi4dKdxXTEOTJyNs/IqVTiQlU68jkuOTPVQ+sPkEuw/Vte8jv0JDkxUda6jxwfdJJW23lU6Ac3ca5NSpQ11zcveD5UqNIkUYixKbanCMG2o77U6Fol6V3DuSExc/apt2b/ZQb3LYj+7on+FYGapP9a5zhiYNiSEENI5GHnNd+zcT2UWmDcC96E2cJRl5SZ0x+b+VDpas/ukIjIVlfnICWpZQieiMmnoAr7Rdmu7v+sUodyL6uWm0o6ptGFzP2twkKZDUT+BThsGmawwbWjqCGv8vbMwmNlUoY+89psl0BSS4X8nqnV/fO9PMP7eWfW5DIo6ei7rCE+bNKSwdm5/bUa5EE10DxpRENSy2qcRVaWOj1FRXRCt9WPWSJ4vShVqpkoFxudtMWlM+dxk1qB4zWeidJ7rxqAchZKhchEqG7wTKRhVvwqFzOfsIqu8dh2W6husND4PYG3zytIedEAwXsiStafG69RflPCuSnt8Soh6CJmBFwioeg/0GCeBqn9lKEvvxjS5oCxdvhUH6vrXQFbi0GRV48pM+GJKMa4Z74H3FrpUoXUYAlWNqyyzsO7mNob1LqPWG/t6gVZ7fOBCLO1jcecLa2N+fZhOlPA9cq5TdYyy7dtt2p0aLKeIa17qedqEUChEnYRpQ0IIIZ2Dkdd8Rw+00u5D8SlEF30F+6NOCSJ8DIQ/ne0xLr0nav6rtoHR9eBjbQwxftE76MR3Cymr52hEVdFy9Tz+XzemrH7J7hboICFOpamIQnQ6LhhDJm59WRqYIgvbPdYp2rLI3esZO7xYTaMS3kMRnx7UqUIBXBRXlhmKwrjHoh2GJQJTRiOqbLk/KdNK023YNGT4/SQ8Rj2HST7nFNFVr4HJdn2wf5QOTLkMSSeheJEKZaGv6l9hei8YwOzaM2Uwzs5m/PFivNXemLB+5Z4LjS8q263HoIQgCwXTPY2eSt5b4+Hchl4A9X7hMuwTIZhk0d+Culmsff3w4lnaTh72sd8WNJUtwk74pclU6jNztyDLxKdoJZpIUqHFS+AF3KYKAVTCVeQoJ+v3p8iCQdPGCZlJpABrAS5Ch6F2HoYDmaVxDr/NXnSUNlSpRlN68TGlH3AciJx1DirBMonPTmO5DB8HXTVKv+zWl3V+VbsS9X5qmZw5MG1ICCGkczDyIk2kffxXtV39Yg1mva23ayNGUVTOQz2uRw9YLvyvWZPZCKSOxFz0lgXGj9QYKzc+qPDtotoGJjfGbKk2UkGDWB1U2l3KKL1Y+icSNbWIGKiWVhlKaENLAXWznGOyNNKIuCy6vZQey1WWBmXh04blpHIYFsbdD502NAVgCuU2bGkJZeKIqiV1GroNBamxYc6soVON+jkL/6aEUZiKtgoVUbWN7XLRWiqKakkZ1vtLcI62foZMNZ4pULyIJ3YIAr7/oXYfppyHonoemtJ92Rp7XvtlYQcg6+cDms7DoPuGPkYPKpZmz0P1ZSnKDi5qvTtvnTZMTltv/PdX2Oncpw7d8wb1PHuMShsWgCBD6b6tc9gLF9G1sP6+GAPxKjI/KLkwkMlMCZZRgqX6GTZSgCYtUnEK0O0jiXP4ZS9+SsgKmx6059YiJaEwarEp4YRIuw1b04bxYy1EalmUwPVC9PEUrjMKihdpoupfrbMu2/10bcuOExP15W7Fzn5Dq/qXUbW0Rv1LfP3KjfsCqmkw3BedFqvw3FXdDW6jq5/5IU5Nq3wpftJLHRnEzXzLFmEr4V6nGHGd+P2hdYSkozDxaqGnnun5PSkGpSjBqiMtmawNGq5+5cUrEDJlh3cCFb0+e3+0QMViFda5JL2sbPM6yg1qW0X43gctpcqyP5NGsKyiqGDAnap5lWq93b+tzkXROiNhzYsQQkjnYORFelPXv1qnTsltU1gJIh33oKw72dpp1XX9Sx2TbNirilNucsHaEm/XQ/1A1p04XMcNhLUxHdHB1qimSBXauQbtPjBAVvg6lfsJWKjiWDXFYr0UTnQoyFDa9FWeuUki9dH1TVE5PLW6NH4gsq5rFVXUZVSdyyirvKtzBVGTaab9EinEZnQVpRSlZdk6SG3KsF+HIQBTR2R9OQzdzZFEJOZTiKmpT1zURTdhp6B4kTSJ+lfP/epJK4G6/mUb6Oa1eOj9U+PBoi+jtrm+UMLPoKs6x1fd7tVll1LNGAybHrRfmv4Yqw3+GFV3ykyw3gqUy6jqmYNt3U3dr0qubApSQgGTqtZlT+jHvfnn9wYNdf9VutRNLBmZMkyh0oZaPIIUohYVuJSeOyZpgW9a4400twXL6seFTRmGtS273DRp2PV9NeANlusU4FR1rh6pQNa5ugHThoQQQjoHIy/SmzrUCKzzKedhMPOyNMwbqUHOsX3eRSDOeWgjCB+FGTFBdOX7IdbpvYThorpO+0RhdKYHLQf2ep0qzBLpRHc+467bBEGSjrWqFKIzcKhIR0rjI0kDfw/ciRKRlzJfQKI0oYqwgsHHOh1YhOsbDkMdGPdwFMZpRHeMPb6IDBqlj8RC5yCaDsN6eVoOQ3d8tC22xQNJazxNGt2C4kVmTq/OG3H9y9noo3RPwm1ol0099keUbR5GYOqxYWKyYMLK6gvaHmOCL1fvdoychzpNpupcxog/pmjWvHRWL6sPKqPOwGLFF1Kl+mwNTMIUondS+udvoJyVUHNzmdK4e2BKU6UDtTBpUar3y6IUYlibQiRYCYFKuA19Os8fE44Zk3BbIUrYVDqxEP+DyJ53KodhWbqaqksZtnXSIHMGpg0JIYR0DkZeZGrUuK/Aedij80bDvGGHNYmBKdR0yPARWdB5I+jSITD2+bMs7Ty0sy3bCMLAz/VVwJ2rmn25fmpn0EiMATP+GBjjIhMbeYUpQvtKTCP6qp+03suO+1IpxEy9BnUpjQhM4C+uh4sQ8YBjFWFlQXrR75PpfobRMakBy0ZiY4Y+PjZv1NFVUaUAbTRtz2G3pQYiG5GqA4s2WaTSho0UoDJclBIe36t/IU0anYLiRfpDWf9mVP9SgtNa/4q/mKJBy0D1RedrY1ACU32BesFSomDEOfK0KIlNNarOHnF6sNpPiVVRPdY9e/09SqcQK01SRnh9ggzeyaiFLCFeWrD8QGQtHCasWSW2NdcjchW2CFEyBdjcr1ELc22f0HQI6nSg3a8sg+Xgc6CFLGq+q4VIn1t6ipz9QSasc3UQpg0JIYR0DkZeZMa4vodAc/BybN7Qo13bzBu6bVRq0LJdtqkno6yDddSlpz7R5g0/H5g/l0sZ6ghrKvOGqdJsymzo74c/pCWFqKIqFzWkx5NZ56FLkfrMbRRFmSjqiaKoIKWo1qsoLCvibdL7GLucMGkYdS4oI0Y1BYraT5svgoa7CD8TvUwaOgVoU35F2dukwelN5gwULzI9RH+LTtH3EHD1r77m/YrThuqLx4pCda5alMoSpv62d+5Clc8L61f1etUFw1hRc5fQX/0rThvqLKaVK4MohRjrmBMl1WUk86Ikma+LAbVgKZFz9niv+a6WFdjg48HIaIoV6roXUNemWgVLiVW/da5YrFT6N6xt2XOXfqaBKG3YqHMlU4g2bdhS59JvAa3xnYbiRabPTOpfOopKTZ2if1XX5o1k096g5uUHZpmiTFvngfbOG7XyuDpZH/UvVNNkBu72VBSmxSpDPWmmqIMspQmHCxi1XlGNk0pY4JWIOFNGr2jJ7RctT3FMqq7VWufSkaS9zqJ6f+0wh0adq2Us17TqXPU+cQuonnUuey7SOVjzIoQQ0jkYeZGZ02ahT9S/Aut8at6vhPOwtWmvTTFBW/JN2joPBIOMJbC9Sx1hqWRfW/0rfOEuJSjSkkLMfcrOndBFSyoK0+fWvRqTVnl/2XGdC6jShSaOolyaLoyo7D426spc2i5MNaaa7zq3oYuG/fl0RGXKaDkajJyscxU+hefciX3Wudx6Pd1JLxhxdRqKFzn16PpVy7xfQed5i005xqkkW/MS+C4cKKsaXGSdr57G+BSVUZb1uv6lO264L3UIytZ2F0BKffScmRmg0p6oxrm5l2p8KktdTpiebH9KnY6LJ5JsjMdKCJEWMitcqW29O2xI2DG+V52rvvbGeK5UnUvC1OK06lzuPkV1rng8lz2GdBqmDQkhhHQORl7kxOjXvFFHUW3zfkkBmFxFWrF1Hs1fyqaOj6oHeZX+c9OlqIBGuwhLtVyv9yk6tc0dh2Bd4ga0b48NGspk4eYdMwijsD6obpvutlGvryOyVregjsJ0mrAlVahNGZluslvaCE1azmcjWx8dmTod2I9JQ0da0zZp1KlFmjTmPhQvcnKYov6lLfFB/ct1uE3UvwovWN4SL8F6mLw+bVlVonRXdl3/sqjMpJiq87sTU3j3on/G2ureU6DsF2o06EuJlRFfFgq2GQRC1hP90lTar9GmaSrnYGyHj1OFkbBVywk7fMs2xPvZ9YGYqVSychQ2u3C01LnU8lTzdNESPzdh2pAQQkjnYORFTi+xecOGQnneNG/kuT9GpZ8k92PHjE071lGXK+gr80YQ0KgHxqBqkeuiMVHhj4+2qrFcfaQQRXz0FRg0aleiDQAyv02Ux6RvRB0TpQl7pw3DtJ/eJzRf6G3NdKB1F8ZNd91ybNKo7402aYQpQHWMdgrWn5WkSaMomiaN+vjYpEHmJtOOvH7605/i2muvxejoKIwx+P73vx9sFxFs27YNo6OjWLhwITZu3IjXX3892GdiYgJbtmzBsmXLsGjRIlx33XV48803T+iFkDMA5w6rXF2uxlDXJAIbs3us/qSyPEtRqnOJTyuV0Xr1Z+y5neOucq1Vg1zVfqX6K+LH9Zdv4eeesl/mWWH/xKfaoj+7PStEPe61b/v2Xn+2C7y9LidAZa/zSbCfFi1bv9INd+3rd9fo7g/C+xffQ/V+GvW+2TShcSLW/l6iKMPt+vMRNd11nyld69KfJbXMlOHcYtrideTIEVxwwQV46KGHktvvv/9+PPDAA3jooYfw8ssvY2RkBFdddRUOHTrk9tm6dSueeuopPPHEE9i1axcOHz6Ma665BkVslSaEEEISGIkbfk3nYGPw1FNP4frrrwdQ/RIaHR3F1q1b8cUvfhFAFWWtXLkSX//61/G5z30O4+PjWL58OR5//HHcdNNNAIC33noLq1evxrPPPourr756yuc9ePAghoeHsREfw4AZnOnlk1OJGzuVwWQGrg9iZmBsOtAY50I0eeb3yTMgz/34qzz3Y77yzJ1b8ixYL3kO5L4hsFseyFxLKMmMW1/mlZlEBrwz0u2XA+LOpdZn/q86h/HjuYxfL1l1Dvuc+hioZYnHdrVlJSW9HE5VEqYKg9Req3kjTBU2XIUqVRikHSfL0MAxWSaX9VguM+nNF4FJQw1MDubvKooqkrI/avVyKRC9PhrX5bcx4jpTmZTjeAFPY3x8HEuWLJn28SfVsLFnzx6MjY1h06ZNbt3Q0BAuvfRSvPjiiwCA3bt34/jx48E+o6OjWLdundsnZmJiAgcPHgz+SIfRKUCX8onSO/UXV5wW0ukk50yrv7Cqwa7ip5K3+xVxWgs+tWi/iOsv42QKsZFO1Gk6nVqL0nuT1Z9JpREnq78s+jNtf+p8WXDu5l/j9bWkOa1wBSnMIHWoUoVFtD64V1Hq0Nbk6vfBvS9BClC9H3GqsJECFC9c9fFuni5bF9NpaNa55gUn1bAxNjYGAFi5cmWwfuXKlfjd737n9lmwYAHOOeecxj72+Jjt27fj3nvvPZmXSk410fiv1tZRbeO/MoQGjl7d59VyW/eNRvf5gMQYMPRYFnhjhqjXqtvNwwdIGWx9Bm6/YCoV1QA46O2rnfuaOPKylKi7WtTbahGzyykjhq+Z9RFtRV00gtZPalvPLhpTzYpcFMp8oepc9qXr8Vx6XdxFg+O55jynxCpvogErItJYF9Nrn3vuuQfj4+Pub+/evSftWgkhhHSPkxp5jYyMAKiiq1WrVrn1+/fvd9HYyMgIjh07hgMHDgTR1/79+7Fhw4bkeYeGhjA0NHQyL5WcLurBy63dN+x+Zdls3tur+0Z1UBAbVQfbKEqFLfHUKfVvtsZPpcAvFDcd1MsIHvtZXUQFEyYILMKalQ6d1DVHFxVcXxRAGH0OfbpSRWNRtBU0342isLB+pZb1enu8chlWzxN1zkh10dDuQaDnVCcuuhKbSlT7JebpCm6TSz8z4prrnNTIa+3atRgZGcGOHTvcumPHjmHnzp1OmC688EIMDg4G++zbtw+vvfZaq3iRjtPri6Rn/cvWNXRNo4z2i+pfde2kUf+yywJv11apRttENrR/J2o/kZ08S9TBwrqQ+tN1LlX/yial+ivEL8d/hfqblKBuFdas0utbrzOu6U2qWlZc1yrC+9PLEh8v2/fXDV8Ihkn444PhE6LShv3UuVjrmldMO/I6fPgwfv3rX7vHe/bswc9//nMsXboU73vf+7B161bcd999OP/883H++efjvvvuw9lnn43PfOYzAIDh4WHceuutuOuuu3Duuedi6dKluPvuu7F+/XpceeWVJ++VEUIImbNMW7x+9rOf4bLLLnOP77zzTgDAzTffjEcffRRf+MIXcPToUdx22204cOAALr74Yvz4xz/G4sWL3TEPPvggBgYGcOONN+Lo0aO44oor8OijjyK3Fmoy93DNUadh3rD5L2XeCCYqSZk31OPAvGGXy1J1Lax3taeLr9mmxfLQy25KCWzw/jnRbuQQnx3U9njV2rDlIhIXqneX5nY9k3HSKh/3Nqwjy+px71ShXa9nRdb7xWlD14eyLH10DATzdgVdNUoBykItN5vuAi0mjXqZKcP5wQmN85otOM6rwxg/5muq8V+uzmWyanxXvezWx+O/ssy3jsr8MXrMF4xx47eq/U04BsyNITPh2K5c7dMY9+WXnXFRjf+y60XPG1bTt3DF6DpXQsCAWIhahMymRJNuw0i8ghpX6feb9ILRGMsV17j0XF1qnJa49XqfuttKos6lx3/RXdhNzqhxXoQQQsjpgI15yemldh9WRGli+ws7y6plG6lkJXpOnWLPG48dsmktwDkMkaNyH9bnilOI6ehHuw1RpwFVSrDeZsRHYRnEXaeUqCMvcacTHYRNJ+KKabgQgdh9WP2ru2iEacKwma/ATPr0YJwqBOCiLttJw/WUBEI3qPjUYtCX0u7n0rpR5BS7C1XaMJiny90DpgrnIxQvMmtIKT3rX/aLqvpibek+r2soiIzsmbfXu5pXARg771jdL14LmDQze9BrjNT1LvtdKT5tVwmXqqDZ66/FSo9j1AOTeyEt23sNXvYd5yVd/1JpQt/EuCVtGKUKq/XKSVjvF7wPqbZP9seEGowsajkYpKwa7zZShak6F5mXMG1ICCGkczDyIqefqHUUgPTgZfuLO8vC1lF68LI9p43I7GDkUnyaTs/IbIwapFy66Kt6nLlMppRhUOQityxhurBRmW4jJaLchaaOvnyk0NZwJo60+skoxoaNpAMxchEGsx0XLQ7DRKoQQHWsSJgqjNODdn08ENmaLKL0IlSEl2oB5c9hXas0acx3KF5k9lD1LykzGF0CU9b5oPuGnjYnh69/JabTcSnEPPNfrnlWpSVRpQ0NfD1N90BE7r8bDfz3Y2DTr49yLwfGT8pYenehMVIJcuAyTMtSY+1U6hW4DiVYZ8pwP+0o9C5A6S9VWErUp1DcfWwMUUiKkhqEXJ87sMRbUdKCpwYju5dBSzypoXiRMwdXqK/qXyY2cKCOwux+Sq8kz5sGDlt7gYqYALePM3JYh35RCVi1zU9xXM3SXJ8WsckiMnNYq72pRMutLkPBMqlBW5qZmDikFrBU5BWYONS4LOkz2tJiZWdOTozZMnELJ222sF1RgEq4dP0qFblxVmTSA9a8CCGEdA5GXmR26VH/SnXfQFl6P18OXyQqimDikmoHHenU9TSocxVlaKPPJNjPugONlH5gM+roK3yiej9fZzOZUY5CWydLR1vOrZiKttoisMSptMuw8XTKmVc179WREqZOFeqmuq7HpEoPJizwccPdYK6teLltIDLrXKQFihc5M4jrX7GBI5VCRLthA0VRrQNCw4aohJ1KH1bnC80bTkyUqPkTpV+GSy/qDKYVtF6vP96eTSNvqIWpYd6Q4LHfT6cDfaqxkSpUKTyfapSG4BjdOSOqc1WrJaxTqXm6JE4vBuLH1k8kDdOGhBBCOgcjL3LmoCx9Uw1grtar8KZ2b7jf5ZnxkUymfqPlmYsSJPdRXP203nlYljC5C51CI4ffrT7GXrdKU2bi3YbuP22vW1+AfTkziDCcOcNGLvH2RIRm3YHOiRilCtsGIus0YpwejE0agE8hpizxRZEeiNx4fYy6iIfiRc48pARM5tJFvVKIQRcOI3AWRJUebOhGpkUl8v3ZbWo8GHK4dkqCEkaMT+uJb+ZbjUGznXmB0OKnffLT+BKeYgbypCClnifaZlRqrmoplXAVioSpwiA1qMd9lU1hc/tFacJgnJd3GLa1fQrqXIQomDYkhBDSORh5kTMPa94w6reVtKQQ7biqLKsHKtcmjbY5wNxGVM+RqTQi4POB0WBmZNZ0kUFyNQBYEIY8dnoUY8KhYGiPglqJj0uRMmFY2owcwfE2nVfvZ2edRpg2bBg0Uk5Cu6yjq6jhbtA9oy0Cjd2FTBeSBBQvcmbScB+29VNS9RUnYEAlYn7Z16w0tbvQde8ow8HM9Xxgpix9CrHuwuFt7eI11hhfd4K0NtU9VcRW+cayJcrChQOOQ4chEqlBJ1wpYdMpwFTD3V6WeIDuQtI3TBsSQgjpHIy8yJlLPw18o/Ffrgei+P1g9LL69V8dCJ8zM+FgZkucQtRpSDvfSb3VtYDK4ihvmkxl1LA0UoXqFPG2VBQTORNND/NFMJar0VjXRsllMoXo3IU9xnNVx0epRkJaoHiRTjGlhb70QgS7nxi4FGQBmFxZ7d12VGlC1w8xawiYe54SQQ9DlyvUglVG4tOvGPlXOs39EdagYvoVtakEyy5HqcKgq4ZuuKsHIseDlgk5AShe5MxH1b8CEwfgvxwzBGPARPw0H4DqtoESYkthWd0FXnfiUMKWjMJcZBXa7etTK5HyX85VNNbnl3W/ifwe7vGwq8ZU9a+E4FkBj00ZjcjLP3aNdVXnjGDCybIMHwfPzToXmT6seRFCCOkcjLxIN2irf6lIrIq2fNrQ1b9KwDkPjbe9A3kVsdm6TJ4H8361phCN8ZGPTge2zdGVckq2pRFPZDxuD6t8Y1tLitGUUQTW2mQ37p5RLxelr1/F0VkcdbHORU4AihfpHu5LLwvqX0EzW20i0DMxK8MHUABiIHXa0MTNfHUKMRYrJT7SZuPX6cTomJ77p5jOl/oUqcIgHZjaJ5na03b6IrS39zk3V3i+KFVIyDRh2pAQQkjnYORFuoVIGPVE/Q9b5wDTKcQgWoPvCJHn6RSijoiyZuRl+hmNHEVVrdHaCRKk/YD2CCsVXcXr9WMVRbmoS7sKZzI3V9t1EtIHFC/SPVT9q9FCaqo5wFIpRJdFLMIUonMW6udoildAL1HSglem1/dNP1/2bQIVb0s4DVPnSo3LqrZ5h2E/Y7mC4/p9LYREMG1ICCGkczDyIt1liv6HfrqUKVKI1sARNPMVP5ZLRwapaEzTjzljJlFb8BxTRCq9OmkAzehKbes5eDgei5Wam0udOzU3F8dzkZMFxYvMDaRsug/tpqlSiMExhTqn7aKhhColaFPRS8hSKck2+nnOHsLkV/Ww1ANTu//i+lVkiU+mHjk3FznJMG1ICCGkczDyIt3GRhFxdCPKzFFKTxdiYOCwgZfJ4Fo6ZSoysdEY9H59kkgLtkZwjdfTHnX1leoD0lFPr5ZRvdBtoKLxXD1NGr2el5BpQPEic4O6/iUqXehSiCZr9kBUx4kbmOy/UI1pCpY7V0A5jVqVWjZKBKcjgDH9ChKQFqWZpPFKSaQeI0t842mEdS5yUqF4kblD1MA3GANm0WPBjGm1houOfESLS1kdp1GNfvvGhOd09HOOtsiolxC1HDPj7u7RTMjJ87CLBjmFsOZFCCGkczDyInOLxADmag6wOqKJBjIDURpRY2tjhXIgGtOc3MTt138UE6Qlg+ufuV2+d+1rGpFPv9FY6pzsokFOExQvMneRVArRqFZGVTPfti/9KUXNEolbX5fWdr5pCCCA3qJ0IiI0DYI6F7tokNME04aEEEI6ByMvMjfRDXxl6hRi8hRtm4qiadqYKVkWRm6W1PlnMkg5wYxNGq3PyS4a5PRD8SJzFz0GTKUQw316OfR6CJsbA3aCIqaEazqC2LcA9TNm62SQmliSkFMI04aEEEI6ByMvMvdpmQNsaprpvLj5b2KX6WNTmk0f48yYhcgneU+ZMiSnEIoXmV+c4Be7FDixjhhJToYCniAnU/AoWuQ0wLQhIYSQzsHIi8wPTmY0IKcoUurXsMHIhhCKFyFnDBQlQvqGaUNCCCGdg+JFCCGkc1C8CCGEdA6KFyGEkM5B8SKEENI5KF6EEEI6B8WLEEJI56B4EUII6RwUL0IIIZ2D4kUIIaRzULwIIYR0DooXIYSQzkHxIoQQ0jkoXoQQQjoHxYsQQkjnoHgRQgjpHCddvLZt2wZjTPA3MjLitosItm3bhtHRUSxcuBAbN27E66+/frIvgxBCyBzmlEReH/jAB7Bv3z739+qrr7pt999/Px544AE89NBDePnllzEyMoKrrroKhw4dOhWXQgghZA5ySsRrYGAAIyMj7m/58uUAqqjrG9/4Br785S/jhhtuwLp16/Dtb38b7777Lr773e+eikshhBAyBzkl4vXGG29gdHQUa9euxac+9Sn85je/AQDs2bMHY2Nj2LRpk9t3aGgIl156KV588cXW801MTODgwYPBHyGEkPnLSReviy++GI899hh+9KMf4Zvf/CbGxsawYcMGvP322xgbGwMArFy5Mjhm5cqVbluK7du3Y3h42P2tXr36ZF82IYSQDnHSxWvz5s34xCc+gfXr1+PKK6/EM888AwD49re/7fYxxgTHiEhjneaee+7B+Pi4+9u7d+/JvmxCCCEd4pRb5RctWoT169fjjTfecK7DOMrav39/IxrTDA0NYcmSJcEfIYSQ+cspF6+JiQn88pe/xKpVq7B27VqMjIxgx44dbvuxY8ewc+dObNiw4VRfCiGEkDnCwMk+4d13341rr70W73vf+7B//3780z/9Ew4ePIibb74Zxhhs3boV9913H84//3ycf/75uO+++3D22WfjM5/5zMm+FEIIIXOUky5eb775Jj796U/jj3/8I5YvX45LLrkEL730EtasWQMA+MIXvoCjR4/itttuw4EDB3DxxRfjxz/+MRYvXnyyL4UQQsgcxYiIzPZFTJeDBw9ieHgYG/ExDJjB2b4cQggh02RSjuMFPI3x8fEZ+RjY25AQQkjnoHgRQgjpHBQvQgghnYPiRQghpHNQvAghhHQOihchhJDOQfEihBDSOShehBBCOgfFixBCSOegeBFCCOkcFC9CCCGdg+JFCCGkc1C8CCGEdA6KFyGEkM5B8SKEENI5KF6EEEI6B8WLEEJI56B4EUII6RwUL0IIIZ2D4kUIIaRzULwIIYR0DooXIYSQzkHxIoQQ0jkoXoQQQjoHxYsQQkjnoHgRQgjpHBQvQgghnYPiRQghpHNQvAghhHQOihchhJDOQfEihBDSOShehBBCOgfFixBCSOegeBFCCOkcFC9CCCGdg+JFCCGkc1C8CCGEdA6KFyGEkM5B8SKEENI5KF6EEEI6B8WLEEJI56B4EUII6RwUL0IIIZ2D4kUIIaRzULwIIYR0DooXIYSQzkHxIoQQ0jkoXoQQQjoHxYsQQkjnoHgRQgjpHBQvQgghnYPiRQghpHNQvAghhHQOihchhJDOQfEihBDSOShehBBCOgfFixBCSOegeBFCCOkcFC9CCCGdg+JFCCGkc1C8CCGEdI5ZFa9//dd/xdq1a3HWWWfhwgsvxH/+53/O5uUQQgjpCLMmXt/73vewdetWfPnLX8Yrr7yCv/mbv8HmzZvx+9//frYuiRBCSEeYNfF64IEHcOutt+Lv/u7v8Bd/8Rf4xje+gdWrV+Phhx+erUsihBDSEWZFvI4dO4bdu3dj06ZNwfpNmzbhxRdfnI1LIoQQ0iEGZuNJ//jHP6IoCqxcuTJYv3LlSoyNjTX2n5iYwMTEhHs8Pj4OAJjEcUBO7bUSQgg5+UziOABAZGZf4rMiXhZjTPBYRBrrAGD79u249957G+t34dlTdm2EEEJOPYcOHcLw8PC0j5sV8Vq2bBnyPG9EWfv3729EYwBwzz334M4773SPy7LE7373O/zlX/4l9u7diyVLlpzya+4aBw8exOrVq3l/WuD9mRreo97w/vRmqvsjIjh06BBGR0dndP5ZEa8FCxbgwgsvxI4dO/Dxj3/crd+xYwc+9rGPNfYfGhrC0NBQsC7LqnLdkiVL+MHpAe9Pb3h/pob3qDe8P73pdX9mEnFZZi1teOedd+Jv//ZvcdFFF+FDH/oQ/u///b/4/e9/j89//vOzdUmEEEI6wqyJ10033YS3334b//iP/4h9+/Zh3bp1ePbZZ7FmzZrZuiRCCCEdYVYNG7fddhtuu+22GR07NDSEr371q410Iqng/ekN78/U8B71hvenN6f6/hiZqU+REEIImSXYmJcQQkjnoHgRQgjpHBQvQgghnYPiRQghpHN0Vrw4Fxiwbds2GGOCv5GREbddRLBt2zaMjo5i4cKF2LhxI15//fVZvOJTz09/+lNce+21GB0dhTEG3//+94Pt/dyTiYkJbNmyBcuWLcOiRYtw3XXX4c033zyNr+LUMdX9ueWWWxqfqUsuuSTYZy7fn+3bt+ODH/wgFi9ejBUrVuD666/Hr371q2Cf+fwZ6uf+nK7PUCfFi3OBeT7wgQ9g37597u/VV1912+6//3488MADeOihh/Dyyy9jZGQEV111FQ4dOjSLV3xqOXLkCC644AI89NBDye393JOtW7fiqaeewhNPPIFdu3bh8OHDuOaaa1AUxel6GaeMqe4PAHz0ox8NPlPPPhv2EJ3L92fnzp24/fbb8dJLL2HHjh2YnJzEpk2bcOTIEbfPfP4M9XN/gNP0GZIO8td//dfy+c9/Plj353/+5/IP//APs3RFs8NXv/pVueCCC5LbyrKUkZER+drXvubWvffeezI8PCz/9m//dpqucHYBIE899ZR73M89eeedd2RwcFCeeOIJt8///M//SJZl8sMf/vC0XfvpIL4/IiI333yzfOxjH2s9Zj7dHxGR/fv3CwDZuXOniPAzFBPfH5HT9xnqXOTFucBC3njjDYyOjmLt2rX41Kc+hd/85jcAgD179mBsbCy4T0NDQ7j00kvn5X0C+rsnu3fvxvHjx4N9RkdHsW7dunlz31544QWsWLEC73//+/HZz34W+/fvd9vm2/2x0y8tXboUAD9DMfH9sZyOz1DnxGu6c4HNZS6++GI89thj+NGPfoRvfvObGBsbw4YNG/D222+7e8H75OnnnoyNjWHBggU455xzWveZy2zevBn//u//jueeew7//M//jJdffhmXX365m09vPt0fEcGdd96JD3/4w1i3bh0AfoY0qfsDnL7P0Ky2hzoR+p0LbC6zefNmt7x+/Xp86EMfwp/92Z/h29/+tiuQ8j41mck9mS/37aabbnLL69atw0UXXYQ1a9bgmWeewQ033NB63Fy8P3fccQd+8YtfYNeuXY1t/Ay135/T9RnqXOQ13bnA5hOLFi3C+vXr8cYbbzjXIe+Tp597MjIygmPHjuHAgQOt+8wnVq1ahTVr1uCNN94AMH/uz5YtW/CDH/wAzz//PM477zy3np+hirb7k+JUfYY6J156LjDNjh07sGHDhlm6qjODiYkJ/PKXv8SqVauwdu1ajIyMBPfp2LFj2Llz57y9T/3ckwsvvBCDg4PBPvv27cNrr702L+/b22+/jb1792LVqlUA5v79ERHccccdePLJJ/Hcc89h7dq1wfb5/hma6v6kOGWfob6tHWcQTzzxhAwODsq3vvUt+e///m/ZunWrLFq0SH7729/O9qWdVu666y554YUX5De/+Y289NJLcs0118jixYvdffja174mw8PD8uSTT8qrr74qn/70p2XVqlVy8ODBWb7yU8ehQ4fklVdekVdeeUUAyAMPPCCvvPKK/O53vxOR/u7J5z//eTnvvPPkJz/5ifzXf/2XXH755XLBBRfI5OTkbL2sk0av+3Po0CG566675MUXX5Q9e/bI888/Lx/60IfkT//0T+fN/fn7v/97GR4elhdeeEH27dvn/t599123z3z+DE11f07nZ6iT4iUi8i//8i+yZs0aWbBggfzVX/1VYNWcL9x0002yatUqGRwclNHRUbnhhhvk9ddfd9vLspSvfvWrMjIyIkNDQ/KRj3xEXn311Vm84lPP888/LwAafzfffLOI9HdPjh49KnfccYcsXbpUFi5cKNdcc438/ve/n4VXc/LpdX/effdd2bRpkyxfvlwGBwflfe97n9x8882N1z6X70/q3gCQRx55xO0znz9DU92f0/kZ4pQohBBCOkfnal6EEEIIxYsQQkjnoHgRQgjpHBQvQgghnYPiRQghpHNQvAghhHQOihchhJDOQfEihBDSOShehBBCOgfFixBCSOegeBFCCOkcFC9CCCGd4/8HEC9nFbRrK3sAAAAASUVORK5CYII=", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAGdCAYAAACl2fynAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAUQNJREFUeJztvX+MXUd5//+ec21vjWWv4jjeH42xrCqoFWulqkMTrJQ4vxz8lfODIGJAqhIpRdDEliwnggaEcKoqhkhN+MNtqg9CCQml5p+YICUCjBKbWlak4AaRpAgZYcBpvHKJzK4d3LV9z3z/OGdmnpkzc+69+8O7Z/f9km587jkzc849e3Oe+zzPe55RWmsNQgghpEFks30BhBBCSK/QeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxrFoti9gMuR5jnfeeQfLly+HUmq2L4cQQkiPaK1x5swZDA8PI8t696MaabzeeecdrFmzZrYvgxBCyBQ5ceIErrzyyp77NdJ4LV++HABwPf4/LMLiWb4aQgghvXIRF3AYL9nnea800niZUOEiLMYiReNFCCGNo6yqO9nUT0+Bxj179uBDH/oQli9fjtWrV+Ouu+7CL3/5S6/NfffdB6WU97ruuuu8NhMTE9ixYwdWrVqFZcuW4Y477sDbb789qQ9ACCFk4dGT8Tp06BAefPBBvPrqqzhw4AAuXryIzZs347333vPaffSjH8XJkyft66WXXvKO79y5E/v378e+fftw+PBhnD17Flu3bkW73Z76JyKEEDLv6Sls+IMf/MB7//TTT2P16tU4evQoPvKRj9j9fX19GBwcjI4xNjaGb37zm3juuedwyy23AAC+/e1vY82aNfjxj3+M2267rdfPQAghZIExpXleY2NjAICVK1d6+w8ePIjVq1fjAx/4AD7zmc/g1KlT9tjRo0dx4cIFbN682e4bHh7GyMgIjhw5MpXLIYQQskCYtGBDa41du3bh+uuvx8jIiN2/ZcsWfOITn8DatWtx/PhxfPnLX8ZNN92Eo0ePoq+vD6Ojo1iyZAkuu+wyb7yBgQGMjo5GzzUxMYGJiQn7fnx8fLKXTQghZB4waeO1fft2/PznP8fhw4e9/du2bbPbIyMjuOaaa7B27Vq8+OKLuPvuu5Pjaa2TqpM9e/bg0UcfneylEkIImWdMKmy4Y8cOfP/738crr7zScXLZ0NAQ1q5di2PHjgEABgcHcf78eZw+fdprd+rUKQwMDETHeOSRRzA2NmZfJ06cmMxlE0IImSf0ZLy01ti+fTuef/55vPzyy1i3bl3HPu+++y5OnDiBoaEhAMCGDRuwePFiHDhwwLY5efIk3nzzTWzcuDE6Rl9fH1asWOG9CCGELFx6Chs++OCD+M53voMXXngBy5cvtzmq/v5+LF26FGfPnsXu3bvx8Y9/HENDQ/jNb36DL37xi1i1ahU+9rGP2bb3338/HnroIVx++eVYuXIlHn74Yaxfv96qDwkhhJA6ejJeTz31FABg06ZN3v6nn34a9913H1qtFt544w08++yz+MMf/oChoSHceOON+O53v+uVAHnyySexaNEi3HPPPTh37hxuvvlmPPPMM2i1WlP/RIQQQuY9SmutZ/siemV8fBz9/f3YhDtZHooQQhrIRX0BB/ECxsbGJpUK4npehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMbRk/Has2cPPvShD2H58uVYvXo17rrrLvzyl7/02mitsXv3bgwPD2Pp0qXYtGkT3nrrLa/NxMQEduzYgVWrVmHZsmW444478Pbbb0/90xBCCFkQ9GS8Dh06hAcffBCvvvoqDhw4gIsXL2Lz5s147733bJvHH38cTzzxBPbu3YvXXnsNg4ODuPXWW3HmzBnbZufOndi/fz/27duHw4cP4+zZs9i6dSva7fb0fTJCCCHzFqW11pPt/L//+79YvXo1Dh06hI985CPQWmN4eBg7d+7EF77wBQCFlzUwMICvfe1r+OxnP4uxsTFcccUVeO6557Bt2zYAwDvvvIM1a9bgpZdewm233dbxvOPj4+jv78cm3IlFavFkL58QQsgscVFfwEG8gLGxMaxYsaLn/lPKeY2NjQEAVq5cCQA4fvw4RkdHsXnzZtumr68PN9xwA44cOQIAOHr0KC5cuOC1GR4exsjIiG0TMjExgfHxce9FCCFk4TJp46W1xq5du3D99ddjZGQEADA6OgoAGBgY8NoODAzYY6Ojo1iyZAkuu+yyZJuQPXv2oL+/377WrFkz2csmhBAyD5i08dq+fTt+/vOf4z/+4z8qx5RS3nutdWVfSF2bRx55BGNjY/Z14sSJyV42IYSQecCkjNeOHTvw/e9/H6+88gquvPJKu39wcBAAKh7UqVOnrDc2ODiI8+fP4/Tp08k2IX19fVixYoX3IoQQsnDpyXhprbF9+3Y8//zzePnll7Fu3Trv+Lp16zA4OIgDBw7YfefPn8ehQ4ewceNGAMCGDRuwePFir83Jkyfx5ptv2jaEEEJIHYt6afzggw/iO9/5Dl544QUsX77celj9/f1YunQplFLYuXMnHnvsMVx11VW46qqr8Nhjj+F973sfPv3pT9u2999/Px566CFcfvnlWLlyJR5++GGsX78et9xyy/R/QkIIIfOOnozXU089BQDYtGmTt//pp5/GfffdBwD4/Oc/j3PnzuGBBx7A6dOnce211+JHP/oRli9fbts/+eSTWLRoEe655x6cO3cON998M5555hm0Wq2pfRpCCCELginN85otOM+LEEKazazO8yKEEEJmAxovQgghjYPGixBCSOOg8SKEENI4aLwIIYQ0DhovQgghjaOneV6ENJYOtTUnP27Df//pfAbGbNzsG9JAaLzI/EcarmkwNiqbIUM4K/iFAXQ+RcOjc3e/acTIDNLwn42EEEIWIvS8yMKg9Lh68pqaHhKcBCqs0NZlWNF6bCqbmVAkIQE0XmT+YsJXoRHqxijNYGiw09p2k2VmKr2V1qxDOFFlhcHSuRb3N2fokMwYC++nJSGEkMZDz4vMTwLvJhourPGupuwdZZf+d2HHK84nH87TdR9HeGUqU77og+INMkPQeJH5RURZ6BkulXlGK2mkJmN8ZkqOP13IJYd6NCYqYvhsmDJTQF7eL53b+61zkf9SigaMTCsMGxJCCGkc9LzI/CQl0gi9rizrzWOaQjhwpoQakimLNlKhxYjXZrwx75yh2lA5j4yQ6YTGi8wflPKMlg0XBobMGhFjuDoYpI5GZyrKxKnK8QOjoDpnvnxCFWGwmnnSGOa5vW8qz11OTKgNVZb7EnqqD8k0wrAhIYSQxkHPi8wPYl6X9GpK70hJT6vc9jyrTl5UJ0/pkpeOCmcV1xCbq9VhUrLnyZX9tdbFPTQhxixzIcRMnEdlkflfZR96YGSK0HiRZpOSxAeGqxIqNNuS0ODJ/R3OG7+2ORDYkMaolbhmz5AIaxYau9IQqTzMc6neQoimDw0YmQI0XmR+oLLuvK2wm1LOOEkZfco4dWvcYue5ROiUMZJ4hkR2Thg7rZ0cPssLA2buqRR5pLwwmDlgFHCQ6WEO/DQkhBBCeoOeF2kuKhHmA6qhQtnHhLiM1yVl9LF6iBHPKulJdeNhTXf1jUDenrw26ZEFIUTnraXChrn7qZtn1vsCSg8sF7msaAgxJqFn/otMHhov0jwiD2cbMkyF8IQhc0ZNVUOFwpBFjUC4L2WIegkTTiak6BmiGtFGJ8NQGh37WQMDZ41aLoxNBmvAAFRDiLFzllU4PAGHgfkvMgkYNiSEENI46HmR5mJEGuW2dygiiXf7A4GGCt6bdilVYugpJesj9uhRdeuBaY0uyvAWpJYyMZ6O8dpink+e27No420VB5z3BaRDiFK8oZ18vvgnp3iDTAkaL9IsRGgvpi5MSeIroUK7nQgVhsavbi7YZPJf06k+rAu5qQ7Gy5AyciasCPj5q0gODIiEEMP8V3QOGPNfpHcYNiSEENI46HmRZiBCe6m1uXoKFZo2daHCLPDQwmtJvQeg60KGMzXnq8ZrUaFnFbaVHpo9Jn7bRkOIgYAD6BhC9CYwm1PL+V+sf0i6hMaLzH3qJgyHE5G7DRWK/tFQYRYaOX/ysmeculEl1n2O6aLmoa/D8KFoq3Ltrk2Okdn/FMfb7WITiRwY0DGE6E9gTuS8uIAl6QKGDQkhhDQOel6kGcREGp3UfHWhwvJ4V6FCpaBbwgOR/4bb4bm9z3AJSkQllzBJe146CBmqtiz3ZPrDKRPb7biAw7QXIUS049ejlPL6e+INqg9JF9B4kblLLM/V5dpclUrxdXJ481D2cmOuj25lvtGS7eR1xq69RNcZrl7jHzXPdpUyXpWclrieXPvGrDTUFSNm3rZa1oABNSHEtvs72PxXrB4iEJfQM/9FaqDxInOTDqWfiiZVgUY0zwWk5fCtlvC8arwtr08WNViegQovfTpl8xmSD3UNlZiz5c5RMXBKA+2qRdStrN4LC3NgQFkGql2OG8l/dVpGxbuu0mLSgJEIzHkRQghpHPS8yNwjtkZXzHOKqQtTea6UHL6XUCFgvS7rZWXV6w0/QzRkOOWfjeWY3YYQRdjQc2TyciTjZQYeWG0IUYQAvRCiFvembhmVRAHfSv1Dqg9JBBovMjcJ81wpAUQ3FeJbrXo5vDBmOiaPDwyZbgXhxbpwoRlL/luDrmkSLZYRnk8njBSUZ+isYVMaWgOqjPShlbmOwmDpVubmiuUmr+XPAyvPAi2EHV1Vok9I6G3+iwIOEoFhQ0IIIY2DnheZO4TqwohgI7kqcs3yJtFQoe3jxtOBFybHsuGuLP7eu/4Sz4vKwmO9y+ZjjldVXaiE+EH53pq5bVpb58p4XNqsqpIr6wWhlbmxRGhPIRPelxwYfghRrtllPCop4Kj7sABM/UOqD0kMGi8yN+igLuxYcDdSIV5JtaBUJQYqQi9UGJnnpbPMPZ+N4fLeC2NY4hmn0E7VGK5U2DBZX1cObh7s5jq0MxCekcvhwnct+AYr09BwuSjPQInQnhLv4eW53H7VyqzwsFKJHm1fNh/Nf4WqSGEMacQWPAwbEkIIaRz0vMicI1wVuauCu5HlTaKhwor4wvfkkqFCE1pUKngP28cLOwqsNzUF1WHFz4hqGPw5Xt68ZHFMQXpkCkDuvC2t7eAamfDYcnexpXekpJso1YqeR2VUkaXXZEOXqmvxRnT1ZbLgofEis0uqikava3OJfsWwgcGLyN6t4eoUKswyFwYslYbmuV2M4T6LV7BXXFolx9V7ystHGr1AVWgMjgZ8haE1FnD3o51DZZm1dLqN4j2Kzrqc3FzsM4MVxiYqo/eUl874VSrRxyYwA4k1wAL1YXltDB0ubBg2JIQQ0jjoeZHZQyoCE+rCahcRKoxMRFYRDysm0vC2uwkVtvx5XtaTUggEH8VmuJ5XrYDDtun48YvuKYdDB/8C9npUrp24Q8PO2dLlvC43z0tWd3Keks6lJxxcrPxsrcyf6FzO+VLtduFFmfNkKr4GWPJTm2sqPWOuvrzgofEic4+ahSU7TkQW77uqWShyXp1ChQCc4TIPUQXvmLUPSm67zxXSrcHy+ojtiiETRqpyzNyCtqtFqHQ5B7jljpntSghRqhCFkarWQDR5Lnlu1XECczFYNf/l103khGVSQONFLj01ea7isOosiTfjhHO5PINXzZl5Ag1T6qmDt6UzaayKbWt0hBdWCDZQtiuNWfiZETFYk8l/ae+fYhiN4jpFzkuZJUky5TU2hk23dbmSsfk8cGkpJLwwZM6AAUA7yH8Jg+fSZFnnOWAAlNJdLWDJ/BdhzosQQkjjoOdFLi2pPFdMXRiTxJsxigGCnJf01oJQoWwnwobdhAqR+f0LT6x4q7Mg/2W8NS8XBv9fwWQqbUjk5GNt/6PEDnfQeGFaeGEKRSTOekRaXKZGbQjRq8RhQoiZk8DLHKCU2RdvRc4LubvYLCGht8OYkGSwgCXzXwsOGi9yaZAGp+su0hDFc14VgxepqlEYm4Rgo4Mww/aXocHMGSatxIrAwsgVuTCI/u5z6fAWTCls6Dpbo2EMU66dAfUqyWu3vlcGKIg8We5fX52Qw8uBib+DhinSK/NXKG6K+UGitZ/LyrTYDkpISel8crFNhhAXGgwbEkIIaRz0vMglp1akUbwRSX+/nRRiKKWsFNt6VKlQYaTChrYTjquhQt2KhAptH3gelpbemnGEMhcS1C0UfWIhwoTXFYo6OsrjUdYc1KKgoZywLMOJSnTUcl+5KSczm31RFaJpqPwKG9ZLhV+Fw5PUt6DKlZiTa4CZ+oddLZ1Cb2uhQeNFZp5e8lzldrT0UzFAei6XDO/VlX2SxibLugoVuj4mz4XyPMp7WMtQoX34B30mM8+rVh4vG2kl5nDBhgeLEGLZP1euinxbl3km08cPQ9aFEJUwzkrO35JlpDwVoqrI54Ey7ybDfmazLCHV2wKWpVFl6HDew7AhIYSQxkHPi8wspdflhQpTpFZFFp5bsuCuEWXEvDIh2NBKekql0jCcgFz2D0OFRR9TYcO9jwkztBB1FJ6X8Kzq5nzZ+xa819FN3wuzIUMp0nBre7mO2j9x6X0BpYcWeF/2OuVE5mDpFBt9DMOJdiDzNyiP5XCh4VgB32IwVFZfbrcRxXr2OdWHCwQaL3LpSeW55HEvz1WXG/ONVVISH5XAo1DHWXm7G0/K4Ss5LhUaJrFt8kIqCBuG7WJGqyZs6BEYMmvAtCqUhJlsZ87pSkBpGSY01yQmM5tjyJyaUQHeRGaFMOdVbopwYlRCL++9WAOsWsAXhQJRVt9oaz+07OW/aKQWGgwbEkIIaRz0vMjMEJaA6lakUfbxw4M1BXfDMF9kPldl5eOWO2fhSbnzuFWRhbcVCDS8ScpKiB9CD02e03uPLsQb/gGVCn9pOE/MeFr2vXTRlGsoPDKltTdXTIYQi1am1JRyE5lzFJ6cCcUiF4INXTv/K6k+lOINIdjoVDrKnCcUb7B01PyHxotMPzF1YaQobbKKhjFcYd1CIFlFw6oLI5L4ykRk00eV/eTEYmOUggnHngoxkMGHuS3bJwgT+sarQwgxQAdWTkmDpaVRkv/GDZaXC2srqCxiwICKETOnUZkqDFQZElRiikAnCX1FfQhU8l/WSJr6h4H6EEC8+oa5ZlbfWBDQeJHpo4tSR17OSiK9MyDIeYk+UvYeCjTqJPEi56U9IymMR0sarIi3BVjDY3NbCZGGX2GjJucVGq+6W6gj2xpuAGOUbIUNcU5hsLz90IUBM4tRyh8ZXo7JXZpu62K3+ayQx1zOKyqhNwbH+0EiPoM3r0uVleRtaQ8//2WQ+a86A6UUDdg8gjkvQgghjaNn4/WTn/wEt99+O4aHh6GUwve+9z3v+H333QellPe67rrrvDYTExPYsWMHVq1ahWXLluGOO+7A22+/PaUPQuYQZZ5LhgzNd8F6S+W23W88LRMKTPZRbrvMc+msekx7bcqcl5HFl+exfVvKLm1iX1npcRmvq/SkCqm8CCsq43G53Jj29gmPzI6ReLVqjpnjNW3s9cWO22PKTaZOfh7lpgTY/u7eGA/V9pF/k6wI6enyb6bl30L+fUx+UtakNKHecr8K/97m+2Qmqsv+lT6Zn2/toaYmaQY9/0Xfe+89XH311di7d2+yzUc/+lGcPHnSvl566SXv+M6dO7F//37s27cPhw8fxtmzZ7F161a0U3M4yNxHJR4YMsSXetCED0DzcOrwqjwAS0Plzd8qH6imnW4FD1TPYImHdUueSz7QY4Yo6BsxXHVGyRjWWsMVtK9v0/31eAYs+Gwm9OqMWDlu+HcQ7VD5oWCuN/P/flnm//BIvOwPG2W+O+5V94PIfaeCH1IyXEkaTc85ry1btmDLli21bfr6+jA4OBg9NjY2hm9+85t47rnncMsttwAAvv3tb2PNmjX48Y9/jNtuu63XSyKEELLAmBFf+uDBg1i9ejU+8IEP4DOf+QxOnTpljx09ehQXLlzA5s2b7b7h4WGMjIzgyJEj0fEmJiYwPj7uvcgcQXhL6Saq8j4M7xQHhLdWCQtVQ4b+L3npZcF5Wbb4LpwsXoYNRQgMCt55PE8l8EyiXlVlv/FUevC2WjWvbr0wG7qU4T0EIU0/hBh6lu6YvL+irbxvmQy5uj7OExOecLd/R+lRCS/KfYkyzwPr/vsqvDjSaKZdbbhlyxZ84hOfwNq1a3H8+HF8+ctfxk033YSjR4+ir68Po6OjWLJkCS677DKv38DAAEZHR6Nj7tmzB48++uh0XyqZZmqrxZuHRbfV4usWlvTW7HLj+jJ4sZ3JPvAqv9sQGeCpCj3lYMu1kaFC2870F/ursvnqGOaYPX8Zokwi528JtaFXBT5cgNIeU+5gDvGztdjvFpCEm7cG0UeuG9YuV6w0n0GIEov7LNbjsn+2wrDJBS3NweQCluX6X171jbrq84Cd/8Xq8/OfaTde27Zts9sjIyO45pprsHbtWrz44ou4++67k/201slfUI888gh27dpl34+Pj2PNmjXTd9GEEEIaxYzP8xoaGsLatWtx7NgxAMDg4CDOnz+P06dPe97XqVOnsHHjxugYfX196Ovrm+lLJb0iwoUdq2jUTUYuBnB9ZBhJvocJ20X2K1WtW+h5bnD9xYRjbz6X9MqER+RVzjAhN9knE2N7Xphrg8wf2x1T/nnk+SUadmKy1sr3woI2sv6uoTK3S3hakF4ZRG1Dzyv074f0kLzeuajQoQG7QnJZ/9BO2ZJ/X7n6cqwCR/mvAmqXTgFKYUgvS6cA4NyvZjIjOS/Ju+++ixMnTmBoaAgAsGHDBixevBgHDhywbU6ePIk333wzabzIHKOLPFeSTnmuaK4j8Qok2lI5WMlzxdSFJq8lJfE2j+Tnuar5pGp+KyWDj6sD3flqpfGJY951dplP83JWsRyZ8q/Nk9B7+00Orarc9O8vvH1eTrLu7yu/K0Ai9xnPf/X2PWb+q8n07HmdPXsWv/rVr+z748eP42c/+xlWrlyJlStXYvfu3fj4xz+OoaEh/OY3v8EXv/hFrFq1Ch/72McAAP39/bj//vvx0EMP4fLLL8fKlSvx8MMPY/369VZ9SAghhNTRs/H66U9/ihtvvNG+N7moe++9F0899RTeeOMNPPvss/jDH/6AoaEh3Hjjjfjud7+L5cuX2z5PPvkkFi1ahHvuuQfnzp3DzTffjGeeeQatVqtyPjK3qRVpFG/c3JugXUWkYQdVlX91pI/9RS+3W7F2IhwXCCf80CVczcJMbIv+3hwpM0YsVJhV98t6hmHpqMr+IHyotIhuiWK8lRJQIgQWroTs7xchRHmiXIv7o7ywnykhZT0oLcKLEfEGMuWN7pWOUuJvLoQYOlPp4r0Q/RPiDVu815aYonhjvqJ0bTGwucn4+Dj6+/uxCXdikVo825ezsKjJc7kmCjA/ROSk0aKTWGdLGDwTGir7eAtIekV2XX+vWrxZl8v2kcpB5ee5rIrQSdmLMdxDPV+kPEMkc15ezioL2nljyf3KPeAD4xXLK4XFepU0WIHx8vbnsOo8lcOp+ITyUO5XuY4cq/ZXuVsPTOXab9fWoo+2dQeLdq4N8rz4FyiMix07t0V+kbsq8Kqdl+1kH5OzEvu16K+L/faxlud2AUvvUZdrmAXKtNlu3qOw0VzUF3AQL2BsbAwrVqzouT8L85LukMKKRG7BE2lIQiPXcXmUIP/hJffltrsumeeS701/LdpVBBrSEInlTTwPRHhCFam8NEqewRPXErYLxqvsD5AF4pPy+NJp0rm5P4EXJr1Pd6QwMJ73J9y9iFdoPdloH+EsSY9XaSuXB4x3LsQc1nsMvK6IdL4YTntCDDuWWbwyKNxb3IM8WrhXZYoCjgYSCSoQQgghcxt6XqQzKTVWhzyX3R/mueywCe/KeFFZ/JiX57KhRnNe5wHU5rnsuF3kucJtGSoMzlMJFYr9XruYd9Lpp6S5BjnJWEroy/1mOJ0r3/uCaOedU7p14QmldxTkv7ycVef8l5TOa4hclBJeu/h+1Oa/lFgcE8G1yNxaF/kvEz70zkPva85D40W6o9v5XOW22x88lFLrdAVCjIox86pqoLodroos+mul4kbJyL/F2lypPJeX4wrG8EKF1sip6n5p8GKCD/kbIfy9IOdsZYiLN+AbtmQIUQ5nDFmnUKE1cuakcpTgYoWuw/tsLUDZP5iIY2ZabCsxB03GIINjsXW/zLhi3a9CnFIeCxevNN+3cPXl0JiROQnDhoQQQhoHPS9STxn2S4k0OvcXnhbge2tm/HDbhAxD7w2+EEPHfpnLX/oR78YU5PX2xzyNQKRR2Y6E/bzVkiNCjqiwI+V5hYhjGnAqQC0chUx4X+I9UHpgyrloqubzdC3ekNctjrlQoYgmKkDDXUMo3vCFHLAD60w5haIS3l4QfraebBuBJxYss1Qj3jDfcYo3mgGNF4kz1TyXNFhCXu+NnSq+a/tUQ4B+2Se37VXSKPukFIbuoau6z3MFbfzQo9lWnlEKc2HR3Fi3vwlkOy3CjsKmyLQTEAkhemG/mhPJgVKGXYT6wvyXfN7Lgr9ATekouS0tXuwHClAWFxbhQNvHGB5t21dKRwHx/FcsXMj815yFYUNCCCGNg54XqUeujNypqedFiXld4t+Kt+Y6Oy/KeGFhuLBs54kyDKa2oPBuYgpDWTw3HjJT/jHEhByIhyQDj6zihUWEHV44sAcvzCoMzWcF3FwuT7Ah+iRCg5U5W6L6hrKaDHHSyH2LFzRWgJk3pzUUxIRw5VaRkGrDWOUNE4pVOcSHkN5esGyKCR0Cbt4X4M/9imFFSRRvNAEaL1IlDPNJahSGYTs5VkwKXVEYeuf3Q4BRhaGKTEpOKQw9Cb0xZNWQYiwcGAufRUOKQZ/QkMXX8BJj94BXdsnsDHNe4hb6ZaJkaA6eKrBilGQ/eVITOiw/RKV0FFCcw4zVUsBF8d5I54FAbSi2zXji/lrpvFQeeguKaRE6hJXOF0OJqRWxyvOxclFyLIYP5xQMGxJCCGkc9LyIIyLSsCHDmNowJdKILVdhqFshOVZ8F74Qw99GdVJybP5USyXCWuH7cPKx8NxCT8l6NCp+TumFZYGAIxVmg+90hCj5w186GsYxKM8lAmuuuZyYLD3RMjSYDPvJUKEJ3xmPKtKnrmwUMnGtyn2gpPJQlwIPLT9sxBuPlY2yYb+2+z5muS/eaIdKxHJsb90vqT5R9L7mEDRepIqUxgfhQy9nFekHoKIc9CppFDsq29FFJjsoDJEF+a9QYSgrzNuHqzRQkQnHdUau3E6GB2se6KkCvJ4hg9gfwbaTz1O4a1F5eUj+6USoTxooFXwee1qNjvdAZ8o3JDXKw7DyhhQSdqU8tI3h/bhQoSpRbCugqLphjtn5yyKGKMctlYeedJ75rzkPjRcpEHmuypyuXippBGNZQ1ZXSSPcznxDVN0uNrVnlFARaXQyKs7zUradHDuZC1MQ15M4pyeh94/JYH0l51VjvGKJLu1uQWGE5DIosoyUGFrDOVTaeFfy2sLcFoBwnldtniy6X1WNpt2uEW8AwsiF3zFpJMX3s932BRsuOea8x7DyRkU6LyqScP7XnIQ5L0IIIY2DntdCJ5LnKvYHUvcUkcUo/eFVuCN97vB9ltgf66PEkvRI545ik5KjYbuUZxG8D3NZKS8q9LBiUnkVpGBS6JboJ/JVJnelPM8pMYaMyKWuO7wHYX/pVUnlYe0k6A54XpQuonzSO7eSeFG0t6I2DIr22rFVuvJG5TqE58bJy3MSGi9S0GmdrlQljcRYSZGGGaPEE2mI47rOkNVVjrfvRVhJhvm864QNZ9n3ESphx/DflCFDvE+lnJJd5BGoe+57uS256HhYH1e53fYhLs4Zhvy0EgYnYYxV5TNLKxngGTVzAWXYMRNv5d9Ki3Zw/QtjHIgmIu3g5biEGEMaLC0MVqzyRtfrfoH5rzkAw4aEEEIaBz2vhYwUVnj7E9L46Bg1CkNzDvlvTFVoqCyfYsJ7SmwnvKgItfL4SLuUPD7Vxn3Y6tjJcGSI9LaE7F0FP+zlZGGtIJYGEfsBT/bubVe8Lf+yokpEef5QNi8c1qQqMfQw5YCh2CYmm/ekiSi+H7JIb+y7k6MaRpTn7DJq6I0b87Y4eXnWofFaiESeTh1Dhp0UhsUgif0quq2zYL/3MBJtU/EBz6iZnFc5dmigAnm8aVPJjaUMHvz9qfPE+lTayo+dw6sQb7fbvvHx+ok5W5XnpoighdfjGShxbRXlYPBZK9vB5/JUiZmuGF7X3p04mRcL8mc272WOCaNky0a1ZQxSbpfGMGawYmt+xYr2JivPc/7XbMOwISGEkMZBz2shE4o0Yiskd+gfrU0o32cRbyzmiSUqakTbSU8mrG7hzSET43ljwW8T8zoQehaIk/DCkh6ZeSumHlmtgfTC8sL7gl3l2bWDho1YKXMsdW1iWwXH7L/K3SMl1ljxRB5hX+ltBcINO1ZswrIZKDVhWRTP1UpVw4hyzldskrIJLdYul1Jez1SVh4rijdmExmuhUZfnimFChqm1ury2yrXxzieMkQrKQEUNWfFPemKy3C+NXmBIOoUAK33iYcewTZ20vJs8l8rhjI8wWDL/Zco5eTmwyHieojALdwTXErkfcpKzee+pEr39wliFhi1xrys/AIQCVEXbafH9MQZJXES41lf5uZW88cGPFSXaRe1UqTxMFu1NVt4IjCJDh5cUhg0JIYQ0DhqvhUT4izRSBioaMixFGtbrUpkLGVrPStYhKj21bkOG5bacZGyVYUa4kSkvZBhSp0B0S5+4WoZRsQYgPATxCvZ7/UNPK9wH917pMsKWO4/LiDTsqw077+tn//Cvbog8eGnXF7kbsxZx/d79iX3myOdOhU/9e6r82o+9Ir9TUrQD4XXH1KpyO5xjKDHfXe87Va+u9b773lAqCLsrhP+PkZmDYcOFhDcztFBMef/z5aZ6AXwDluc2pKJyWJVWsdCfbQSXoNFOudVqFeOaIrlaxLUCSbNdUBBleCgXT327RlXagKUoQm5CwSZyR5VAj8wr6fh+JY7psE04YdjsF0ZDZWJoWX9Qu8tEG/jLrz7ghghqJXrFhMX+WgNmolzBh5afp+5zJwWC4rqVTqsNuyLX7qYa2bs9j/ZveGpbKgHz4GJMfkpr107n8bW87NDmeoLcXtiHYcNLCo3XQkPLbL+cJJQjmvcy8X9btDRzD6dW2DZi/Mz/0PIhJNq5/SKfkSugBWfMNLxSQFo8UJX3cA0MU+zBW/vgFtp0YdjCNlqeR44prkfL/aGnk7ldXloquGbVhqukIbyeynIrwvjFPqPdF7kfMUOmIuMU+7XXRiXuQex+uD5mDPFhxQ+X4l4HFyUXnYTYNn1yfywERi5lfNz4umiTS8MWl8ob4+cZLs1812zAsCEhhJDGQc9rIaPdWkUqUy6kkmeujl5dDF/nTm6cAe4nckt4UbkLHQJF+DDqbWkbWlTGu/LkbiKUZKTYGt6S94WEXIv3Yjxv1rAYM/DWJEr8mE/WHAw8Femt6Ujo0L7NXDvpUXnXIkKN4RpiMlSYJPB8wmP2X+ERyWMV7wp+H+ft+IeV3C//vvI+5oGHZ0N42h9HeGV2nKBdNGQovSjbLg+up8v4pgwvSiiRn1VovBYiYdkdoJr/sk11UW1AyIij+S8ZQ8y1b/wiVb+L/totGClDPUqVRso8aFQ1RBl8FhvaK5/24TPXhjqlVSn7eG0T4S/veGiwUuEzkbKT1d414Bss2Tco92S3gzyXJ0GPnT9yPTFD5u5b+rNWtmVfBOHAxPPctHGGrSbMJsOEcrwgHKi80F18W8vclkTkuWzIUPzgihXmLbrFxmLIcDZg2JAQQkjjoOe1kKkTb9hwYJfSvlz71S68c4h/pUcl9yslQn7a8wptGBHwln+vhBYFFUVgFx6EHwpz4UQdaaOjmnE3tnd+4x3FLlWGChU8RWFl+EjYMCa+qHhWQTvPcwzaqcj9KfbrivjCXRii9zcZdrT9jNesxXUGYULRDqGnJb9bYTgx5g2lvLA6UuFFijRmHRovUhDmv+QhEzoEulMeSuPXUl7FAq+Uj1JOeSgNVvkAc6kxJaTyIpSklAvFmQoUpk+uhaRP+TJ8e51AqDCMEaoAK//GDFb4cBfXVon0hWrBDoS5sYqBShixtCGTRqGqEPQ+h+2v09crDJEzKuVLFCH2/lbyxw1c/5Q8vjZkKOXxKeNjFIaACxmmFIZm6F7yZGTGYdiQEEJI46DntdCJiDeK/XIichehw9DbssMHATYt4mdSyCGvxRyX3tY0VS6ITlgWzkFSsJEIkyXDk3KsUqwhD3nzvMztaIkDMWJjByTDeal2OtIudQ+C/rGJyZWQ32SIeWFyf+x9sB31nqZDYZi6FnLJofEiBSL/pfOsUnnDLUFfozxsVcdCnrnQIVCGHcv/8eUk5VDurJzR0yK0VWxXJywrU6nCPmy7nLBsQmAw+TRx0kCVWFwbyrifuR7/PEq2EyQnJgfvawnCgOF+VRb29SrTy8K+ud+uGEf7hihHEDbUYtu00Wnjngg7GqWhErmp2MRkGSa02zKfpSPbXtgxd/sAVHJTocLQjCUVhqlKGsxzzSkYNiSEENI46HmRKp54wy8bVRFvBP0A+DUPVSp0aA8WW8K7C+fzIBeTpLxtiHCiBtpSoyHVi4CcsOzVTJSeU+6crVrPQnhVcv5WysMzxXit9yqGM+/dgG5cSd0kYxVsd1IYhrUMpVdW97lDFaGc21U7Mdnsb0vPCb5HFZvbZb4HwttSoaowsq3DY8lJxs5bqy0DRZHGnITGizgi+a9i8nKNdD6mPMxkkd7cr9jhhR21kHy7vIlWOl5xo7ItZPNQXsUNJUtXiHBiVDZv32soa9jE/YiFIEUIzZpmqSiUBg6l/RfP7VgIUVIrMRcdpSHyqszDDxtWDbAzPPWhPjm29seqMWxuuzxPWwfXmpDHh1J5ERKs1DD0jJ9QCoYGx/6oChSGIcLIcTLy3IfGi1SROauwWG83+S8Ir8eWjRIV5+V5YqKM8Bc34HtLkWrzyqzMG/MuxAKEWjnZvPXiUrkt2weVB3V0blgujH7m57+kwYKoKq/MGEDF20oR5rtkvkqWXfKMhZf/0tU2EUNUyZ8FRskTaUhD1HbGRgX5Knfvxd83R3w79j2I9Zf7AVTKQNkmzshFDV4M5rnmLMx5EUIIaRz0vEg9Zf5LZQlPLEau/TW/Ws4TK9acEjmz0NtC8Svfhg4BIAvVhnK7rGWYF55bZdIy4OVXVKaqeamYFyVChYX34IcGKwWBy/6ht2Vvmeijys8E0RXoIlQYwYYKy3NW3scmBQceZkVtGIQUiz662j/yuVU7aBee34wt7nVUYQhY7yr03iyxnFUe/NthrS4P5rkaBY0XiRPJfxX7Xf6rYohSlTdkGDKsOG/Gl/OvTBWOiLAjJd5QZRtpmPzQnjmPtqsuW/FGOO/L9pHb4pyAi1l4zziX//JyXHlgsMz78qPZeV69hA1l2E8YCC/nFRhTaZRU0EeGAOP5ryAXFog0lDBKdWt2yTxXUqQBsR2MEctzJUOGQLySRrntGTzO52ocDBsSQghpHPS8SD2l15Sqe9i5f+5Ch4Av+AjqHPqS6DJ0WPYxP7NSykNb2Ldura9yLBupVH6FjVB5GBVy1ITMZFhMF75gQVb1tuRl2wOT8Lw8UUXphYUCjmI7ItJIfJ54GDXdR+XSIxLt8vi9Kd4nFIaBdzWpZU+SxXSdtxarwMHJyM2Cxot0h80HiPxXUgIfqbyhxRMtt5bIhX7C4rsQfUJVYkx5iBwqy2yfpPIwmNeFinBdWhBrYkQTXRlDtPKk/7q8NrPf5rwyiM/jbKRnyOoIjIIXGhTGDCK0V5cLkzL4iirRk92XbdrBsTCkGKoF7X6IChgQITy5HYQJIzmwYjtQFJqbq8vtWCWNyn10xop5rubBsCEhhJDGQc+LdEa7cJy/v0vxhrc+GJxAQYvKG5GQkLdcijjmqQ3NeVBOrLIeTWfloam8UZn3VZzUD5MFwgw7YEyYYcjK68xVVbAhQ4UR5WEU4UnGJh9XBRu+5xQTchiPqluRhtn2RSw6eq88j9d4TfZ6nOcUCjsqk5IjxwpxTiSMKD0ugCKNeQyNF+kOEfarW/cLKB9MLbHAVyidV267tvKGfGjJsaTaMHMPQCmd70p5WD5cnYTdGSV/YrL2+igbboRtD8APDcobk2lnwMxBM5y4HXXKw4rCUBifUHnolW0KZPThfhv+C4yZvD+ACxXae1AxcnFDFBbf9XNbENtxA5VSGHp5LimHj0nju12ni0arUTBsSAghpHHQ8yK9oV1cKBRveHUPZagmOu8L5TiibFQo3jBzuNo5dCtz7YSHpyHFH2LelygJBUT9JEd53Z7IQ05MDkaorGQso0/iMv39QsARCDOigo0az6tSEsrbH4QKI6HGmHelpBfUyQszocZYGai225ZtrKclRBpKztOS3lpb7E+JNAB4Ig2DmdcVmwNWHjd9qS5sNjReZEqYB4As3hvLf1n1YZj/kpU37KDVcKGX/6qpvKGN0Wqh+/yXfd9F/iti1HRp/DxjVX4mb7+1OPD7CINVMWSJ0Ke3HYbzcnnMXUM0T5b7RqpjngsQRkWMV5fnAgqjFqgPYyrCSvHdVJ4rD/vLcGLnPFe08C5pFDRepHek59RNuSiJzH+h5X49e3PBcn/RypaqPtDK04fVNqSivav8Vy4tBFyFjJSE3nhdFQFHaZCkR2SOGqW9l+cyxlj54o3UnK/wuiMGxraRwoyYFxZK4ANPLpTEA/ArdLT9Chud8lxem27yXAZ77ry6L6waPxljFPPcSGNgzosQQkjjoOdFJk9N/quiIozVPVRa/HzKoc2Ckd748H4ZF+tniZCk+GWv27AptNSilUV/uPNkQKyShsx/eaE5E8vLRB/rBAgPTXpOpUOmgn32HoiKH7EZCRbhbdXJ5sO6heFnqCxv4nliQThQtvO243kubypC6Gm14xOYkefx4rtB2FBrF/bzPLdY/cKEKpF5rvkDjReZNmrzX6L6hpXR1+W/elm0smgElWX2eeSVkcqL4wCgW1k1/2XateXDTIYTFbQxitYgCGPmFel1IUS3OGd5LGq85LZGJVwoEcajTjYfqyRfqYiRkMpDw4pdQqPmxBdI5rlU2zdElRJQtsKGM0qdFpnslOcq3uc1fZjnmo8wbEgIIaRx0PMiU6ML8YbW2qkPWy0xmRlp8Uanuofm3KF4Q9Y6LH+b6VYg3jDDhteZqUr1jeI6RTiwVBc6QYY4JmsWQnpRyvew3NBFPWKh9k96XuJjd5LNu/civJj77cI1u2IiDeVJ2+VYuhRtuGNRkUZbelqoF2lEvas86BMXaaQmH0ehSGPeQONFpode5n91mf+qKx0FxPNflXJRQGnQiv0KmTe1LCxsK42UPWWmKjkvXx4vjsVCg9CAUr50PqUq7JKKwtBsSxWgpzD0jVqYC5OhQimJd3O5hFELclFKtEOeB4pCM27eW56r3O4pz1VuM8+1MGDYkBBCSOOg50WmHy3rH+Zx9SHQcfJytE87B1quvw0NGvFGbN2vSiWOctNcbvmvaok3WgFCwGE9La2BthNzVEKIxsvMEXhaIiYYzueSzkAXgg3zMcy+cMJyrPpGuOyJHwIEpFrQm89l+wvxRRkylGHDVOWMyjpdUqQhRRZeVY1AKVg3GTnZR4g0uNTJvKRnz+snP/kJbr/9dgwPD0Mphe9973veca01du/ejeHhYSxduhSbNm3CW2+95bWZmJjAjh07sGrVKixbtgx33HEH3n777Sl9EDIHkCGfZJPyQROusWSLo+ZuOxcPL9PHPJzMfvnK/T6qXYaztHYPX/ngLUNkVgYujhkjoDRsfqfoY8JkIj9kwo85/P453EvXHJOvds0r1ScxLhLn9Pa3TRhQ7C/vW3jNaJf3r2xr7rV/zP09zL13f7O8+jeTBkoHbULDI78T3nenCBV2zH0xZDiv6Nl4vffee7j66quxd+/e6PHHH38cTzzxBPbu3YvXXnsNg4ODuPXWW3HmzBnbZufOndi/fz/27duHw4cP4+zZs9i6dSva7XZ0TEIIIUSidE9SnaCzUti/fz/uuusuAMWv6uHhYezcuRNf+MIXABRe1sDAAL72ta/hs5/9LMbGxnDFFVfgueeew7Zt2wAA77zzDtasWYOXXnoJt912W8fzjo+Po7+/H5twJxapxZO9fDKTmNm2KiuWTjFKxExBmWNZ5tplmdvfahXhN/HeHlPKzRNTyoUQlSqK92ZZ5ZjOMvczLcugW+6cWsEV/c2UDQ8ic6FCnSk3j7lVtjEFeJWYAyZqFmqhIjRiDS1ChXJbomtmKavY/6ph2FBLMUc4t0uEBo1AoxRVREOFXmhQ9M+L4rlO6JJHBRsqzyvhQK/orrfdttv2kdRu+yKNdrtepFFeG0UazeCivoCDeAFjY2NYsWJFz/2nVbBx/PhxjI6OYvPmzXZfX18fbrjhBhw5cgQAcPToUVy4cMFrMzw8jJGREdsmZGJiAuPj496LzHFsWEgssV6GfWyIxwsRidCPDR268JE7ln7ZPEoQnlLmgVrKtW04Mc9dtQjtQmA2PBgLIRojUI5nQ29tBGG6MKzXRZgwB7K2Tr5SfdJhxGCf1y/4nMHnk8et4TJhQl0auPJe++HBwmipMERo1IUd/oY2Z2XDxe6744WcxbbfJ/i+0XDNW6ZVsDE6OgoAGBgY8PYPDAzgt7/9rW2zZMkSXHbZZZU2pn/Inj178Oijj07npZLZxMzPCavPy+OiSgeU265WrC8fTuXilV71DfPLvpW56vNhCSkxH0xONSt0FE5g4aTyKAoF288QyOiN49TSbs6ZAiDnfSk3DSBcdNISKcwbbSb3ax14Xs5bM8YVEHk6CMMlljdxVeEDb8uOZQyHO68Rt3hFetvOeFRWRfaEGIkVjmXuMyTsQxYcMyKVV0HoQ2td2RdS1+aRRx7B2NiYfZ04cWLarpUQQkjzmFbPa3BwEEDhXQ0NDdn9p06dst7Y4OAgzp8/j9OnT3ve16lTp7Bx48bouH19fejr65vOSyWXCq0BBNL51DIq4eTlzLpRsD/zTfUN09784MnL/xgPToaLcm37e/UPjYReLmApJjB7y5q0ZLUN+J6R/dElvCsoIJPemegj3CVRyKMnVMzZKMN+8aVPpLcF5xEJVaEdw/PCnEcVLiwZlcTnwlOSnpZUENoxZD7M7BehPhFmtmOnwoChNJ7hwnnPtHpe69atw+DgIA4cOGD3nT9/HocOHbKGacOGDVi8eLHX5uTJk3jzzTeTxos0HJPH6jH/ZfMYXq5E5saCl2in2v54dRJ6I/OWx2SYS8rH5VwqG26z7QBfcu7nl9J5LvdKy+Fr2nn5L5Gnk7kt7dp7n9MYLnPdclqAlmOJ+9zWXk6rIomX+cZ27v4W8m8S5rlkjkrI4XvJc9FwLSx69rzOnj2LX/3qV/b98ePH8bOf/QwrV67E+9//fuzcuROPPfYYrrrqKlx11VV47LHH8L73vQ+f/vSnAQD9/f24//778dBDD+Hyyy/HypUr8fDDD2P9+vW45ZZbpu+TEUIImbf0bLx++tOf4sYbb7Tvd+3aBQC499578cwzz+Dzn/88zp07hwceeACnT5/Gtddeix/96EdYvny57fPkk09i0aJFuOeee3Du3DncfPPNeOaZZ9AyEmgyr6ldOgVw6391s3SKkVi3ZMHCAiXaeSICW6GjFGyUoUYp5jBScndt4gPIeIUSJTIywMQGtQJUZs6pSi2HuB6Jim7WE3EuvHW9ACuDL7YD8UXu2lhVIYJQofG4gLgowxNsyBBgXhnL87aAQPYu9oeCjQ6SeMB9n8jCYkrzvGYLzvNqMMrN+Qrnf7km5VwuUzpKKXdcZWJulyv6q7LMnzfW8ueQuflbmT83LJwDZvvLOV9iWwXzvDLl+iiXG9NizlddEd66eV3dUJn7pf1tV0wXwhAJFaKQyxfHgmrxdo5VmNfyc17SeNm5XJ7hEeE/ANrLjeWiv9vvhZPL8ZLGi9XiG8ecmudFCCGEXApYmJdcWkr1IYDJLZ3Sgvh1LeeCle27WX1Zzv+SRXuDEKIs3CvnZSlxCQruoM6U825Ch0p4WFJdqDyFYg/EwoZmnwyd2lCf66PEnC1Xq1B6YmJbhgrF2lxKekRyPpfY9lSEUnQB+OHBQBlaUReGRXcBK9Iw2/S4Fh40XuTSI4xPWH0eEPmvbqrPG0No1v+KLWApJPSqnbtyUGX+y1sDTErqTQml4mBBOZFZmTJQbQ0v8me2MwV1sQx/mQFMiFPeiqlFDauSeZH/kcc8qbxRB5o24QRkz7D5ocJiLKEeBPxQoWfkZC4rd8pCwDc4MgRYkdAzVEjiMGxICCGkcdDzInOPXHcuHQWUocZyf7sN3WrFV1/O4U8MLr0E3co89WFdCNFTIcoyUqIAr2oj/nPQTqSuhhQn44VFJygDlVCiFHNIRaEnyijFG0lVYRgqtG2C8GAsVBjs1+ExGfaTisFQoEE1IYlA40Vmj1T+K1aBI5b/kvUPy3qDzpC0nXwe8PJc/vnhHo6Ze+Dr3FX1kCFEtEojVVoQ3cq8qhumHqISOTcFXVUVRgzVpCOINUbLhU59RaELIeZezisZKpSGKKhT6OW5KjmvIF8V1i0st708V/JzMs9FHDReZHapyX91XH05R7WEVGreWOV8Iv9lH5hBzssU1g29sEzb92jnNpdVmQ8mLq0QZiTM01SC9+GzXn6+8NlujYrzPl11eNMn8LY8Q+bnr7w8V8wLiwk0ZM4rlMSb/t2sikzDteBhzosQQkjjoOdF5hw6152rbyAIISrtqQ8rBXylpyRCh0qEGotxRc6r9JQ8FaJSkfflWGVrAF6Y0G3qcIc5VfUehFHGbhyNmgnLfijRLTgZnXAcqgrLPtI7qlTPiKkNtYaWk4896Xyi4G7dwpKECGi8yNxA5L8qOS+ZywpCiJUSUrZrDq2FND21Blju1vPSmWgHN5Yn5Mh05L35DO6SizyX+AyZJ81wW4lQYi/5r0qVjYRc3p1UhvZCIxUasqqBsYt+Rss7uf224K4cSwtLXRcqjME8FxEwbEgIIaRx0PMic4dAvFGrPrRdjDcAX32YASiL9noSek/NllnvCwAUMj+EaFwnMZEZWkG34L+XYUFxqd7iqnlcsKECqWCdXH5yocOE8jAQWHhS+YRyUIVeWFg9o2xnw37tdlWwIUUaHT9LTpEGSULjReYmukP1DSCe/4rNAZNhRyBuwMr9qrQeWkrq2zlkwV/VFkYmC1YAtzk4V21ey8UqxWljTEoun1CXe+HElFEzisJU/iomh7chwESosGYuV7d5rqJLEGYkRMCwISGEkMZBz4vMPbQfYrPqQ6Dj/K/oBGa4Ehk6gx9ClDUQZchP1kCUmHqI9j0QUxgic55OUby3RmE4mWVR6kJosWOxkKFsn5izFZ3LZbwus5aaDBV2mojcSaQRUxcyXEgi0HiRuUkH9SEQDyFGJzDLek7S+AHFA9goFqVSHkEZKYmYmFw0LFWNXpjON0i1pZ+ySTyca6JpFeMEdMiFBWHDulChbBMNFebVsW2bmlBhDKoLSQ0MGxJCCGkc9LzI3CVVOipV+zAUcJjgnlzbKywhBbjwV6vVXQhRa9+VioX9Mn9fNHQoLmnaSHkqdaG4YLs2VAhYFWE0VBh6Yfb8RqGY+LAsAUV6hMaLNAP7cMvi+a+YYbCTkQMJPfz+Ln/VFnUQs2gIEUpBZ8p/sMbOXat5n+IiXnWkHvh1YUOERXsTocIwxxULFaYK7lYuJ5HnorqQdAmNF5n7BAIOtz8ioe9UwDeowlHJgYUyevGPuRZVqd1UvTadzYLxqvFUVCcRRMwTkyWhYnL4lLdVU3A36q318DkIMTDnRQghpHHQ8yLNIJr/8itYdFXANxZCDHNgQMQDkxOMO3teKo/vr+szaXqVzcf2xxR/qcoZ4fImdWtzherCyHmY5yKTgcaLNAshoa+UkKop4FsfQkzI6JPigs7GC0DVyHWiW4PWywM+KUPvYi5VsnJG3lOosNK/7MOFJclUYNiQEEJI46DnRZpLWP9QSOi11t2HEOtk9IBT2QFV7yjLakJzc8Dz6tQn9C6DdrVy+G5Dhd6ArFtIpgcaL9I8zAMxVkLK5rNUvE8shJiS0Zs8lxxLB4V4pWELrqmWLBL0mGroLBXm7GL8qKQ9CPNVclypUGH02rQfKuzyughJwbAhIYSQxkHPizSXWP3DWA1DKd6QIoK6ZVRk2QvpXCnhkRkCzyxGZRJ16LFNI12tlWWozP+qmXtVt7yJN2ZCpBGehx4XmQI0XmR+oM0ClsECkEBVQi+7aZ2oRJ+QxuvS6Hihw8Q1yRxcaPAuFd3klbouKVVTOaNTwd1wYUlCpgjDhoQQQhoHPS/SbALxhhEEeOrDmvlfyDLrKfjLqAiCsGF54vj1xLy0uUIqhCep84zC/l3O5yr+6WJeGSE9QONF5gcm/1UaF28By3ICMxAJIQpjprUG2tWHqp+vihikrIsQYopYhfxemGoYroNBS+bPulmbK1QXMs9FphEaLzJ/CAQcUS8MEJ5WabjakTyWHFa+iUncIwYvRrTyfc/Wrp6exBoxOsnti5N4bWvncsl9NFxkGmHOixBCSOOg50XmJzpYtFJI6E2Yz6vCERLzsIB6iXuHCcpzyu/o1QuK3CfP40qtx0V1IZkhaLzI/EIKOLwFLEvDItYAQ6aSYbakUeuWlPGbTabwmWrDkYEEPrqwJEOGZJqZg/+HEUIIIfXQ8yLzk2D1ZSfeqE5ijnaH7n1JE0lNeDEu3Jg8UxZpdEM3MntQEk8uHTReZP4iFrCMSug7kc9MYGLWqm1Mhi5zVhVlIUDDRWYUhg0JIYQ0DnpeZGGgq/O/OhMP/amphBPnON3fm+QA03MhhHSAxovMf2T+axoerjpaLqqhzISxYbiQXAIa/n8eIYSQhQg9L7IwmClvYK4V3yVkgUDPixBCSOOg8SKEENI4aLwIIYQ0DhovQgghjYPGixBCSOOg8SKEENI4aLwIIYQ0DhovQgghjYPGixBCSOOg8SKEENI4aLwIIYQ0DhovQgghjYPGixBCSOOg8SKEENI4aLwIIYQ0DhovQgghjYPGixBCSOOYduO1e/duKKW81+DgoD2utcbu3bsxPDyMpUuXYtOmTXjrrbem+zIIIYTMY2bE8/rgBz+IkydP2tcbb7xhjz3++ON44oknsHfvXrz22msYHBzErbfeijNnzszEpRBCCJmHzIjxWrRoEQYHB+3riiuuAFB4XV//+tfxpS99CXfffTdGRkbwrW99C3/84x/xne98ZyYuhRBCyDxkRozXsWPHMDw8jHXr1uGTn/wkfv3rXwMAjh8/jtHRUWzevNm27evrww033IAjR44kx5uYmMD4+Lj3IoQQsnCZduN17bXX4tlnn8UPf/hDfOMb38Do6Cg2btyId999F6OjowCAgYEBr8/AwIA9FmPPnj3o7++3rzVr1kz3ZRNCCGkQ0268tmzZgo9//ONYv349brnlFrz44osAgG9961u2jVLK66O1ruyTPPLIIxgbG7OvEydOTPdlE0IIaRAzLpVftmwZ1q9fj2PHjlnVYehlnTp1quKNSfr6+rBixQrvRQghZOEy48ZrYmICv/jFLzA0NIR169ZhcHAQBw4csMfPnz+PQ4cOYePGjTN9KYQQQuYJi6Z7wIcffhi333473v/+9+PUqVP4p3/6J4yPj+Pee++FUgo7d+7EY489hquuugpXXXUVHnvsMbzvfe/Dpz/96em+FEIIIfOUaTdeb7/9Nj71qU/h97//Pa644gpcd911ePXVV7F27VoAwOc//3mcO3cODzzwAE6fPo1rr70WP/rRj7B8+fLpvhRCCCHzFKW11rN9Eb0yPj6O/v5+bMKdWKQWz/blEEII6ZGL+gIO4gWMjY1NSsfA2oaEEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIaB40XIYSQxkHjRQghpHHQeBFCCGkcNF6EEEIax6war3/913/FunXr8Cd/8ifYsGED/vM//3M2L4cQQkhDmDXj9d3vfhc7d+7El770Jbz++uv4m7/5G2zZsgW/+93vZuuSCCGENIRZM15PPPEE7r//fvzd3/0d/uIv/gJf//rXsWbNGjz11FOzdUmEEEIawqwYr/Pnz+Po0aPYvHmzt3/z5s04cuTIbFwSIYSQBrFoNk76+9//Hu12GwMDA97+gYEBjI6OVtpPTExgYmLCvh8bGwMAXMQFQM/stRJCCJl+LuICAEDryT3EZ8V4GZRS3nutdWUfAOzZswePPvpoZf9hvDRj10YIIWTmOXPmDPr7+3vuNyvGa9WqVWi1WhUv69SpUxVvDAAeeeQR7Nq1y77P8xy//e1v8Zd/+Zc4ceIEVqxYMePX3DTGx8exZs0a3p8EvD+d4T2qh/ennk73R2uNM2fOYHh4eFLjz4rxWrJkCTZs2IADBw7gYx/7mN1/4MAB3HnnnZX2fX196Ovr8/ZlWZGuW7FiBb84NfD+1MP70xneo3p4f+qpuz+T8bgMsxY23LVrF/72b/8W11xzDT784Q/j//2//4ff/e53+NznPjdbl0QIIaQhzJrx2rZtG95991384z/+I06ePImRkRG89NJLWLt27WxdEiGEkIYwq4KNBx54AA888MCk+vb19eErX/lKJZxICnh/6uH96QzvUT28P/XM9P1RerI6RUIIIWSWYGFeQgghjYPGixBCSOOg8SKEENI4aLwIIYQ0jsYaL64FBuzevRtKKe81ODhoj2utsXv3bgwPD2Pp0qXYtGkT3nrrrVm84pnnJz/5CW6//XYMDw9DKYXvfe973vFu7snExAR27NiBVatWYdmyZbjjjjvw9ttvX8JPMXN0uj/33Xdf5Tt13XXXeW3m8/3Zs2cPPvShD2H58uVYvXo17rrrLvzyl7/02izk71A39+dSfYcaaby4Fpjjgx/8IE6ePGlfb7zxhj32+OOP44knnsDevXvx2muvYXBwELfeeivOnDkzi1c8s7z33nu4+uqrsXfv3ujxbu7Jzp07sX//fuzbtw+HDx/G2bNnsXXrVrTb7Uv1MWaMTvcHAD760Y9636mXXvJriM7n+3Po0CE8+OCDePXVV3HgwAFcvHgRmzdvxnvvvWfbLOTvUDf3B7hE3yHdQP76r/9af+5zn/P2/fmf/7n+h3/4h1m6otnhK1/5ir766qujx/I814ODg/qrX/2q3fd///d/ur+/X//bv/3bJbrC2QWA3r9/v33fzT35wx/+oBcvXqz37dtn2/zP//yPzrJM/+AHP7hk134pCO+P1lrfe++9+s4770z2WUj3R2utT506pQHoQ4cOaa35HQoJ74/Wl+471DjPi2uB+Rw7dgzDw8NYt24dPvnJT+LXv/41AOD48eMYHR317lNfXx9uuOGGBXmfgO7uydGjR3HhwgWvzfDwMEZGRhbMfTt48CBWr16ND3zgA/jMZz6DU6dO2WML7f6Y5ZdWrlwJgN+hkPD+GC7Fd6hxxqvXtcDmM9deey2effZZ/PCHP8Q3vvENjI6OYuPGjXj33XftveB9cnRzT0ZHR7FkyRJcdtllyTbzmS1btuDf//3f8fLLL+Of//mf8dprr+Gmm26y6+ktpPujtcauXbtw/fXXY2RkBAC/Q5LY/QEu3XdoVstDTYVu1wKbz2zZssVur1+/Hh/+8IfxZ3/2Z/jWt75lE6S8T1Umc08Wyn3btm2b3R4ZGcE111yDtWvX4sUXX8Tdd9+d7Dcf78/27dvx85//HIcPH64c43cofX8u1XeocZ5Xr2uBLSSWLVuG9evX49ixY1Z1yPvk6OaeDA4O4vz58zh9+nSyzUJiaGgIa9euxbFjxwAsnPuzY8cOfP/738crr7yCK6+80u7nd6ggdX9izNR3qHHGS64FJjlw4AA2btw4S1c1N5iYmMAvfvELDA0NYd26dRgcHPTu0/nz53Ho0KEFe5+6uScbNmzA4sWLvTYnT57Em2++uSDv27vvvosTJ05gaGgIwPy/P1prbN++Hc8//zxefvllrFu3zju+0L9Dne5PjBn7DnUt7ZhD7Nu3Ty9evFh/85vf1P/93/+td+7cqZctW6Z/85vfzPalXVIeeughffDgQf3rX/9av/rqq3rr1q16+fLl9j589atf1f39/fr555/Xb7zxhv7Upz6lh4aG9Pj4+Cxf+cxx5swZ/frrr+vXX39dA9BPPPGEfv311/Vvf/tbrXV39+Rzn/ucvvLKK/WPf/xj/V//9V/6pptu0ldffbW+ePHibH2saaPu/pw5c0Y/9NBD+siRI/r48eP6lVde0R/+8If1n/7pny6Y+/P3f//3ur+/Xx88eFCfPHnSvv74xz/aNgv5O9Tp/lzK71AjjZfWWv/Lv/yLXrt2rV6yZIn+q7/6K0+quVDYtm2bHhoa0osXL9bDw8P67rvv1m+99ZY9nue5/spXvqIHBwd1X1+f/shHPqLfeOONWbzimeeVV17RACqve++9V2vd3T05d+6c3r59u165cqVeunSp3rp1q/7d7343C59m+qm7P3/84x/15s2b9RVXXKEXL16s3//+9+t777238tnn8/2J3RsA+umnn7ZtFvJ3qNP9uZTfIS6JQgghpHE0LudFCCGE0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMZB40UIIaRx0HgRQghpHDRehBBCGgeNFyGEkMbx/wMvhDr6/2DLcwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] @@ -263,11 +309,26 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "id": "028ad375", "metadata": {}, - "outputs": [], - "source": [] + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Including gas 0 transition at 345.7959899 GHz\n", + "65536\n" + ] + } + ], + "source": [ + "g = Gas()\n", + "g.set_properties_from_lambda('co.dat')\n", + "\n", + "image = model.make_image(npix=256, pixel_size=0.2*u.arcsec, channels=np.linspace(-5., 5., 10)*u.km/u.s, rest_frequency=g.nu[2], \n", + " incl=45.*u.degree, pa=45.*u.degree, distance=1.*u.pc, include_dust=False, device='cpu')" + ] } ], "metadata": { diff --git a/pinballrt/camera.py b/pinballrt/camera.py index 4b9575f..fd9fad1 100644 --- a/pinballrt/camera.py +++ b/pinballrt/camera.py @@ -52,6 +52,7 @@ def emit_rays(self, x, y, nu, nx, ny, pixel_size): ray_list = PhotonList() ray_list.position = wp.array(position, dtype=wp.vec3) ray_list.direction = wp.array(direction, dtype=wp.vec3) + ray_list.direction_frame = wp.array(direction, dtype=wp.vec3) ray_list.indices = wp.zeros(xflat.shape+(3,), dtype=int) ray_list.intensity = wp.array2d(intensity, dtype=float) ray_list.tau_intensity = wp.array2d(tau_intensity, dtype=float) diff --git a/pinballrt/gas.py b/pinballrt/gas.py new file mode 100644 index 0000000..fdc8563 --- /dev/null +++ b/pinballrt/gas.py @@ -0,0 +1,101 @@ +import astropy.constants as const +import astropy.units as u +import urllib +import requests +import numpy +import os + +class Gas: + + def set_properties_from_lambda(self, filename): + if not os.path.exists(filename): + if os.path.exists(os.environ["HOME"]+"/.pinballrt/data/gas/"+filename): + filename = os.environ["HOME"]+"/.pinballrt/data/gas/"+filename + else: + web_data_location = 'https://home.strw.leidenuniv.nl/~moldata/datafiles/'+filename + response = requests.get(web_data_location) + if response.status_code == 200: + if not os.path.exists(os.environ["HOME"]+"/.pinballrt/data/gas"): + os.makedirs(os.environ["HOME"]+"/.pinballrt/data/gas") + urllib.request.urlretrieve(web_data_location, + os.environ["HOME"]+"/.pinballrt/data/gas/"+filename) + filename = os.environ["HOME"]+"/.pinballrt/data/gas/"+filename + else: + print(web_data_location+' does not exist') + return + + f = open(filename) + + for i in range(3): + f.readline() + + self.mass = float(f.readline()) + + f.readline() + nlev = int(f.readline()) + f.readline() + + self.J = numpy.empty(nlev, dtype=" ksi[i]].min() == cum_lum)] - cell_coords = wp.array2d(np.array(cell_coords)[:,:,0], dtype=int) + t1 = time.time() + cell_coords = wp.array2d(np.zeros((nphotons_per_source, 3), dtype=int), dtype=int) + wp.launch(kernel=self.sample_from_cum_lum, + dim=(nphotons_per_source,), + inputs=[self.grid, wp.array3d(cum_lum, dtype=float), cell_coords, np.random.randint(0, 100000)]) + t2 = time.time() + timing["cell coords time"] = t2 - t1 + t1 = time.time() new_position = wp.array(np.zeros((nphotons_per_source,3)), dtype=wp.vec3) wp.launch(kernel=self.random_location_in_cell, dim=(nphotons_per_source,), inputs=[new_position, cell_coords, self.grid, np.random.randint(0, 100000)]) + t2 = time.time() + timing["random loc time"] = t2 - t1 photon_list.position = wp.array(np.concatenate((photon_list.position.numpy(), new_position.numpy()), axis=0), dtype=wp.vec3) + t1 = time.time() new_direction = wp.array(np.zeros((nphotons_per_source, 3)), dtype=wp.vec3) wp.launch(kernel=self.random_direction, dim=(nphotons_per_source,), inputs=[new_direction, torch.arange(nphotons_per_source, dtype=torch.int32, device=wp.device_to_torch(wp.get_device())), np.random.randint(0, 100000)]) photon_list.direction = wp.array(np.concatenate((photon_list.direction.numpy(), new_direction.numpy()), axis=0), dtype=wp.vec3) + photon_list.direction_frame = wp.array(np.concatenate((photon_list.direction_frame.numpy(), new_direction.numpy()), axis=0), dtype=wp.vec3) + t2 = time.time() + timing["random dir time"] = t2 - t1 + t1 = time.time() photon_list.frequency = wp.array(np.concatenate([photon_list.frequency.numpy(), np.repeat((const.c / wavelength).to(u.GHz).value, nphotons_per_source)]), dtype=float) photon_list.energy = wp.array(np.concatenate([photon_list.energy.numpy(), np.repeat(self.total_lum/nphotons_per_source, nphotons_per_source).astype(np.float32)]), dtype=float) + t2 = time.time() + timing["set freq/energy time"] = t2 - t1 return photon_list @@ -146,7 +234,7 @@ def tau_distance_scattering(grid: GridStruct, ix, iy, iz = photon_list.indices[ip][0], photon_list.indices[ip][1], photon_list.indices[ip][2] - photon_list.alpha[ip] = grid.density[ix,iy,iz] * photon_list.ksca[ip] + photon_list.alpha[ip] = grid.dust_density[ix,iy,iz] * photon_list.ksca[ip] distances[ip] = photon_list.tau[ip] / photon_list.alpha[ip] @wp.kernel @@ -253,7 +341,7 @@ def photon_cell_properties(photon_list: PhotonList, ix, iy, iz = photon_list.indices[ip][0], photon_list.indices[ip][1], photon_list.indices[ip][2] photon_list.temperature[ip] = grid.temperature[ix, iy, iz] - photon_list.density[ip] = grid.density[ix, iy, iz] + photon_list.density[ip] = grid.dust_density[ix, iy, iz] photon_list.amax[ip] = grid.amax[ix, iy, iz] photon_list.p[ip] = grid.p[ix, iy, iz] @@ -377,9 +465,9 @@ def update_grid(self, timing={}): temperature = ((total_energy*u.L_sun).cgs.value / (4*const.sigma_sb.cgs.value*\ planck_mean_opacity*\ - self.mass.cgs.value))**0.25 + self.dust_mass.cgs.value))**0.25 - temperature[np.logical_or(temperature < 0.1, self.mass.value == 0)] = 0.1 + temperature[np.logical_or(temperature < 0.1, self.dust_mass.value == 0)] = 0.1 if (np.abs(old_temperature - temperature) / old_temperature).max() < 1.0e-2: converged = True @@ -395,7 +483,7 @@ def initialize_luminosity_array(self, wavelength): for i in range(self.shape[0]): for j in range(self.shape[1]): for k in range(self.shape[2]): - self.luminosity[i,j,k] = (4*np.pi*u.steradian*self.grid.density.numpy()[i,j,k]*self.volume.cpu().numpy()[i,j,k]*self.dust.interpolate_kabs(nu)*self.distance_unit**2*models.BlackBody(temperature=self.grid.temperature.numpy()[i,j,k]*u.K)(nu)).to(u.au**2 * u.Jy).value + self.luminosity[i,j,k] = (4*np.pi*u.steradian*self.grid.dust_density.numpy()[i,j,k]*self.volume.cpu().numpy()[i,j,k]*self.dust.interpolate_kabs(nu)*self.distance_unit**2*models.BlackBody(temperature=self.grid.temperature.numpy()[i,j,k]*u.K)(nu)).to(u.au**2 * u.Jy).value self.total_lum = self.luminosity.sum() @@ -793,6 +881,81 @@ def check_pixel_too_large(photon_list: PhotonList, photon_list.pixel_too_large[ir] = pixel_size > cell_size[ix,iy,iz] + def select_lines(self, lam): + # Convert wavelengths to frequencies and find range + nu = (const.c / lam).to(u.GHz) + max_nu = nu.max() + min_nu = nu.min() + + line_nu = [] + alpha_line = [] + inv_gamma = [] + + # Check which lines fall in range for each gas + for igas, gas in enumerate(self.gases): + # Calculate maximum line width + max_v = self.grid.velocity.numpy().max() * const.c + a_thermal = 3*np.sqrt(2 * const.k_B * self.grid.temperature.numpy().max()*u.K / \ + (gas.mass * const.m_p)) # 3-sigma width + a_microturb = self.grid.microturbulence.numpy().max() * u.km / u.s + + max_v += a_thermal + a_microturb + + # Check each transition + for iline in range(gas.nu.size): + max_frequency = gas.nu[iline] * (1 + max_v / const.c) + min_frequency = gas.nu[iline] * (1 - max_v / const.c) + + if ((min_nu < min_frequency and min_frequency < max_nu) or + (min_nu < max_frequency and max_frequency < max_nu) or + (min_frequency < min_nu and max_nu < max_frequency)): + + print(f"Including gas {igas} transition at {gas.nu[iline]}") + + alpha_line_tmp, inv_gamma_tmp = self.calculate_level_populations(igas, iline) + + line_nu.append(gas.nu[iline].to(u.GHz).value) + alpha_line.append(alpha_line_tmp) + inv_gamma.append(inv_gamma_tmp) + + with wp.ScopedDevice(self.device): + self.grid.alpha_line = wp.array4d(np.array(alpha_line), dtype=float) + self.grid.inv_gamma = wp.array4d(np.array(inv_gamma), dtype=float) + self.grid.nlines = len(line_nu) + self.grid.line_nu = wp.array1d(np.array(line_nu), dtype=float) + + def calculate_level_populations(self, igas, iline): + gas = self.gases[igas] + + # get level indices (convert to 0-based) + level_up = gas.J_u[iline] - 1 + level_low = gas.J_l[iline] - 1 + + level_populations_upper = np.zeros(self.shape) + level_populations_lower = np.zeros(self.shape) + alpha_line = np.zeros(self.shape) + self.inv_gamma_thermal_np = np.zeros(self.shape) + + # partition function: either callable or array-like + Q = gas.partition_function(self.grid.temperature.numpy()) + + # level populations (LTE Boltzmann with partition function) + level_populations_upper = gas.g[level_up] * np.exp(-const.h * const.c * gas.E[level_up] / (const.k_B * self.grid.temperature.numpy()*u.K)) / Q + level_populations_lower = gas.g[level_low] * np.exp(-const.h * const.c * gas.E[level_low] / (const.k_B * self.grid.temperature.numpy()*u.K)) / Q + + # thermal / microturbulent broadening + a_thermal = np.sqrt(2.0 * const.k_B * self.grid.temperature.numpy()*u.K / (gas.mass * const.m_p)) + a_microturb = self.grid.microturbulence.numpy() * u.km / u.s + inv_gamma = (1.0 / (gas.nu[iline] / const.c * np.sqrt(a_thermal**2 + a_microturb**2))).decompose().to(1./u.GHz) + + number_density = (self.grid.gas_density.numpy()*u.g / u.cm**3 * self.gas_abundances[igas] / (gas.mass * const.m_p)).decompose() + + # Alpha for this line in this cell + # alpha = c^2/(8*pi*nu^2) * A * n * (n_l * g_u/g_l - n_u) * inv_gamma / sqrt(pi) + alpha_line = (const.c**2 / (8.0 * np.pi * gas.nu[iline]**2) * gas.A_ul[iline] * number_density * (level_populations_lower * (gas.g[level_up] / gas.g[level_low]) - level_populations_upper) * (inv_gamma / np.sqrt(np.pi))).decompose().to(1./self.distance_unit) + + return alpha_line, inv_gamma + @wp.kernel def add_intensity(ray_list: PhotonList, s: wp.array(dtype=float), @@ -807,23 +970,39 @@ def add_intensity(ray_list: PhotonList, tau_cell = 0. intensity_abs = 0. + intensity_sca = 0. alpha_ext = 0. alpha_sca = 0. - tau_cell = tau_cell + s[ir]*ray_list.kext[iray,inu]*grid.density[ix,iy,iz] - alpha_ext = alpha_ext + ray_list.kext[iray,inu]*grid.density[ix,iy,iz] - alpha_sca = alpha_sca + ray_list.kext[iray,inu]*ray_list.ray_albedo[iray,inu]*grid.density[ix,iy,iz] - intensity_abs = intensity_abs + ray_list.kext[iray,inu] * (1. - ray_list.ray_albedo[iray,inu]) * \ - grid.density[ix,iy,iz] * planck_function(ray_list.frequency[inu], grid.temperature[ix,iy,iz]) + if grid.include_dust: + tau_cell = tau_cell + s[ir]*ray_list.kext[iray,inu]*grid.dust_density[ix,iy,iz] + alpha_ext = alpha_ext + ray_list.kext[iray,inu]*grid.dust_density[ix,iy,iz] + alpha_sca = alpha_sca + ray_list.kext[iray,inu]*ray_list.ray_albedo[iray,inu]*grid.dust_density[ix,iy,iz] + intensity_abs = intensity_abs + ray_list.kext[iray,inu] * (1. - ray_list.ray_albedo[iray,inu]) * \ + grid.dust_density[ix,iy,iz] * planck_function(ray_list.frequency[inu], grid.temperature[ix,iy,iz]) + + albedo_total = alpha_sca / alpha_ext + + if alpha_ext > 0.: + intensity_abs = intensity_abs * (1.0 - wp.exp(-tau_cell)) / alpha_ext + + intensity_sca = (1.0 - wp.exp(-tau_cell)) * albedo_total * scattering[inu,ix,iy,iz] + + intensity_line = float(0.) + if grid.include_gas: + vector_velocity = wp.vec3(grid.velocity[0,ix,iy,iz], grid.velocity[1,ix,iy,iz], grid.velocity[2,ix,iy,iz]) + for itrans in range(grid.nlines): + # Need to add doppler shift here based on cell velocity along ray direction + profile = wp.exp(-(ray_list.frequency[inu] * (1.0 + wp.dot(ray_list.direction[ir],vector_velocity)) - grid.line_nu[itrans])**2. * (grid.inv_gamma[itrans,ix,iy,iz]**2.)) - albedo_total = alpha_sca / alpha_ext + alpha_this_line = grid.alpha_line[itrans,ix,iy,iz] * profile - if alpha_ext > 0.: - intensity_abs = intensity_abs * (1.0 - wp.exp(-tau_cell)) / alpha_ext + tau_cell = tau_cell + s[ir] * alpha_this_line + alpha_ext = alpha_ext + alpha_this_line - intensity_sca = (1.0 - wp.exp(-tau_cell)) * albedo_total * scattering[inu,ix,iy,iz] + intensity_line = intensity_line + alpha_this_line * planck_function(ray_list.frequency[inu], grid.temperature[ix,iy,iz]) - intensity_cell = intensity_abs + intensity_cell = intensity_abs + intensity_sca + intensity_line ray_list.intensity[ir,inu] = (ray_list.intensity[ir,inu] + intensity_cell * wp.exp(-ray_list.tau_intensity[ir,inu])) * (1. - wp.float32(ray_list.pixel_too_large[ir])) ray_list.tau_intensity[ir,inu] = ray_list.tau_intensity[ir,inu] + tau_cell @@ -838,7 +1017,7 @@ def reduce_source_intensity(ray_list: PhotonList, ix, iy, iz = ray_list.indices[ir][0], ray_list.indices[ir][1], ray_list.indices[ir][2] - ray_list.intensity[ir, inu] = ray_list.intensity[ir, inu] * wp.exp(-s[ir] * ray_list.kext[ir, inu] * grid.density[ix, iy, iz]) + ray_list.intensity[ir, inu] = ray_list.intensity[ir, inu] * wp.exp(-s[ir] * ray_list.kext[ir, inu] * grid.dust_density[ix, iy, iz]) def propagate_rays(self, ray_list: PhotonList, frequency, pixel_size): with wp.ScopedDevice(self.device): @@ -1015,9 +1194,9 @@ def random_location_in_cell(position: wp.array(dtype=wp.vec3), rng = wp.rand_init(seed, ip) - position[ip][0] = grid.w1[ix] + (grid.w1[ix+1] - grid.w1[ix])*wp.randf(rng) - position[ip][1] = grid.w2[iy] + (grid.w2[iy+1] - grid.w2[iy])*wp.randf(rng) - position[ip][2] = grid.w3[iz] + (grid.w3[iz+1] - grid.w3[iz])*wp.randf(rng) + position[ip][0] = grid.w1[ix] + (grid.w1[ix+1] - grid.w1[ix])*max(min(wp.randf(rng), 0.999), 0.001) + position[ip][1] = grid.w2[iy] + (grid.w2[iy+1] - grid.w2[iy])*max(min(wp.randf(rng), 0.999), 0.001) + position[ip][2] = grid.w3[iz] + (grid.w3[iz+1] - grid.w3[iz])*max(min(wp.randf(rng), 0.999), 0.001) @wp.kernel def next_wall_distance(photon_list: PhotonList, @@ -1087,7 +1266,7 @@ def minimum_wall_distance(photon_list: PhotonList, if sz2 < s: s = sz2 - if s * photon_list.kabs[ip] * grid.density[ix, iy, iz] < 10.**log10_tau_min: + if s * photon_list.kabs[ip] * grid.dust_density[ix, iy, iz] < 10.**log10_tau_min: s = 0. max_tau_distance = 10.**log10_tau_max / photon_list.alpha[ip] @@ -1235,6 +1414,8 @@ def photon_loc(photon_list: PhotonList, i3 += 1 photon_list.indices[ip][2] = i3 + photon_list.direction_frame[ip] = photon_list.direction[ip] + class UniformSphericalGrid(Grid): def __init__(self, ncells=9, dr=1.0*u.au, mirror=True, device="cpu"): """ @@ -1339,9 +1520,9 @@ def random_location_in_cell(position: wp.array(dtype=wp.vec3), rng = wp.rand_init(seed, ip) - r = grid.w1[ix] + wp.randf(rng) * (grid.w1[ix+1] - grid.w1[ix]) - theta = grid.w2[iy] + wp.randf(rng) * (grid.w2[iy+1] - grid.w2[iy]) - phi = grid.w3[iz] + wp.randf(rng) * (grid.w3[iz+1] - grid.w3[iz]) + r = grid.w1[ix] + max(min(wp.randf(rng), 0.999), 0.001) * (grid.w1[ix+1] - grid.w1[ix]) + theta = grid.w2[iy] + max(min(wp.randf(rng), 0.999), 0.001) * (grid.w2[iy+1] - grid.w2[iy]) + phi = grid.w3[iz] + max(min(wp.randf(rng), 0.999), 0.001) * (grid.w3[iz+1] - grid.w3[iz]) position[ip][0] = r * wp.sin(theta) * wp.cos(phi) position[ip][1] = r * wp.sin(theta) * wp.sin(phi) @@ -1487,7 +1668,7 @@ def minimum_wall_distance(photon_list: PhotonList, if sp < s: s = sp - if s * photon_list.kabs[ip] * grid.density[iw1, iw2, iw3] < log10_tau_min: + if s * photon_list.kabs[ip] * grid.dust_density[iw1, iw2, iw3] < log10_tau_min: s = 0. max_tau_distance = log10_tau_max / photon_list.alpha[ip] @@ -1701,6 +1882,14 @@ def photon_loc(photon_list: PhotonList, photon_list.position[ip][1] = photon_list.radius[ip] * photon_list.sin_theta[ip] * photon_list.sin_phi[ip] photon_list.position[ip][2] = photon_list.radius[ip] * photon_list.cos_theta[ip] + # Also, since we have all of the components calculated, update the direction frame + + xhat = wp.vec3(photon_list.sin_theta[ip] * photon_list.cos_phi[ip], photon_list.cos_theta[ip] * photon_list.cos_phi[ip], -photon_list.sin_phi[ip]) + yhat = wp.vec3(photon_list.sin_theta[ip] * photon_list.sin_phi[ip], photon_list.cos_theta[ip] * photon_list.sin_phi[ip], photon_list.cos_phi[ip]) + zhat = wp.vec3(photon_list.cos_theta[ip], -photon_list.sin_theta[ip], 0.) + + photon_list.direction_frame[ip] = photon_list.direction[ip][0] * xhat + photon_list.direction[ip][1] * yhat + photon_list.direction[ip][2] * zhat + class LogUniformSphericalGrid(UniformSphericalGrid): def __init__(self, ncells=9, rmin=0.1*u.au, rmax=4.5*u.au, mirror=True, device="cpu"): """ @@ -1812,9 +2001,9 @@ def random_location_in_cell(position: wp.array(dtype=wp.vec3), if ix == 0: r = 10.**(wp.randf(rng) * grid.logw1[ix]) else: - r = 10.**(grid.logw1[ix-1] + wp.randf(rng) * (grid.logw1[ix] - grid.logw1[ix-1])) - theta = grid.w2[iy] + wp.randf(rng) * (grid.w2[iy+1] - grid.w2[iy]) - phi = grid.w3[iz] + wp.randf(rng) * (grid.w3[iz+1] - grid.w3[iz]) + r = 10.**(grid.logw1[ix-1] + max(min(wp.randf(rng), 0.999), 0.001) * (grid.logw1[ix] - grid.logw1[ix-1])) + theta = grid.w2[iy] + max(min(wp.randf(rng), 0.999), 0.001) * (grid.w2[iy+1] - grid.w2[iy]) + phi = grid.w3[iz] + max(min(wp.randf(rng), 0.999), 0.001) * (grid.w3[iz+1] - grid.w3[iz]) position[ip][0] = r * wp.sin(theta) * wp.cos(phi) position[ip][1] = r * wp.sin(theta) * wp.sin(phi) @@ -1972,3 +2161,11 @@ def photon_loc(photon_list: PhotonList, photon_list.position[ip][0] = photon_list.radius[ip] * photon_list.sin_theta[ip] * photon_list.cos_phi[ip] photon_list.position[ip][1] = photon_list.radius[ip] * photon_list.sin_theta[ip] * photon_list.sin_phi[ip] photon_list.position[ip][2] = photon_list.radius[ip] * photon_list.cos_theta[ip] + + # Also, since we have all of the components calculated, update the direction frame + + xhat = wp.vec3(photon_list.sin_theta[ip] * photon_list.cos_phi[ip], photon_list.cos_theta[ip] * photon_list.cos_phi[ip], -photon_list.sin_phi[ip]) + yhat = wp.vec3(photon_list.sin_theta[ip] * photon_list.sin_phi[ip], photon_list.cos_theta[ip] * photon_list.sin_phi[ip], photon_list.cos_phi[ip]) + zhat = wp.vec3(photon_list.cos_theta[ip], -photon_list.sin_theta[ip], 0.) + + photon_list.direction_frame[ip] = photon_list.direction[ip][0] * xhat + photon_list.direction[ip][1] * yhat + photon_list.direction[ip][2] * zhat diff --git a/pinballrt/model.py b/pinballrt/model.py index b315227..a704135 100644 --- a/pinballrt/model.py +++ b/pinballrt/model.py @@ -2,6 +2,7 @@ import astropy.units as u from .grids import Grid from .dust import load, Dust +from .gas import Gas from .camera import Camera from schwimmbad import SerialPool import xarray as xr @@ -50,7 +51,8 @@ def __init__(self, grid: Grid, grid_kwargs={}, ncores = 1, pool = SerialPool()): self.ncores = ncores self.pool = pool - def add_density(self, density: u.Quantity, dust): + def set_physical_properties(self, density=None, dust=None, amax=None, p=None, gases=None, abundances=None, + velocity=None, microturbulence=None): """ Add density to the grid. @@ -63,7 +65,9 @@ def add_density(self, density: u.Quantity, dust): """ for device in self.grid_list: for grid in self.grid_list[device]: - grid.add_density(density, load(dust) if isinstance(dust, str) else dust) + grid.set_physical_properties(density=density, dust=load(dust) if isinstance(dust, str) else dust, + amax=amax, p=p, gases=gases, abundances=abundances, + velocity=velocity, microturbulence=microturbulence) def add_star(self, star): """ @@ -191,7 +195,8 @@ def scattering_mc(self, nphotons, wavelengths, device="cpu", return_timing=False if return_timing: return timing - def make_image(self, npix=100, pixel_size=None, lam=np.array([1.])*u.micron, incl=0, pa=0, distance=1*u.pc, nphotons=100000, device="cpu"): + def make_image(self, npix=100, pixel_size=None, channels=None, rest_frequency=None, incl=0, pa=0, distance=1*u.pc, include_dust=True, + include_gas=True, include_sources=True, nphotons=100000, device="cpu"): """ Create an image from the dust distribution. @@ -222,9 +227,34 @@ def make_image(self, npix=100, pixel_size=None, lam=np.array([1.])*u.micron, inc if pixel_size is None: pixel_size = ((1.25*self.grid.grid_size()*self.grid.distance_unit / distance).decompose()*u.radian).to(u.arcsec) / npix + self.grid.grid.include_dust = include_dust + self.grid.grid.include_gas = include_gas + + # Check whether spectral is wavelength or frequency + + if channels.unit.is_equivalent(u.micron): + lam = channels.to(u.micron) + nu = (const.c / channels).to(u.GHz) + elif channels.unit.is_equivalent(u.GHz): + nu = channels.to(u.GHz) + lam = (const.c / nu).to(u.micron) + elif channels.unit.is_equivalent(u.km / u.s): + if rest_frequency is None: + raise ValueError("rest_frequency must be provided when channels are in velocity units.") + nu = (rest_frequency * (1 - channels / const.c)).to(u.GHz) + lam = (const.c / nu).to(u.micron) + else: + raise ValueError("Either lam or nu must be provided.") + + # Check which lines from the gas should be included + + if include_gas: + self.grid.select_lines(lam) + # First, run a scattering simulation to get the scattering phase function - self.scattering_mc(nphotons, lam, device=device) + if include_dust: + self.scattering_mc(nphotons, lam, device=device) # Now set up the image proper. @@ -252,8 +282,9 @@ def make_image(self, npix=100, pixel_size=None, lam=np.array([1.])*u.micron, inc zip(self.camera_list[device], np.array_split(new_x, self.ncores), np.array_split(new_y, self.ncores))))).sum(axis=0) * (u.Jy / u.steradian) * \ image.pixel_size**2 - source_intensity = np.array(list(self.pool.map(lambda camera: camera.raytrace_sources(image.x, image.y, nx, ny, image.nu, distance, - nrays=int(1000/self.ncores)).numpy(), self.camera_list[device]))).mean(axis=0) * u.Jy + if include_sources: + source_intensity = np.array(list(self.pool.map(lambda camera: camera.raytrace_sources(image.x, image.y, nx, ny, image.nu, distance, + nrays=int(1000/self.ncores)).numpy(), self.camera_list[device]))).mean(axis=0) * u.Jy intensity += source_intensity diff --git a/pinballrt/photons.py b/pinballrt/photons.py index 34f97ca..37afd55 100644 --- a/pinballrt/photons.py +++ b/pinballrt/photons.py @@ -40,4 +40,6 @@ class PhotonList: cos_theta: wp.array(dtype=float) phi: wp.array(dtype=float) sin_phi: wp.array(dtype=float) - cos_phi: wp.array(dtype=float) \ No newline at end of file + cos_phi: wp.array(dtype=float) + + direction_frame: wp.array(dtype=wp.vec3) \ No newline at end of file diff --git a/pinballrt/sources.py b/pinballrt/sources.py index 90fb03f..84605fb 100644 --- a/pinballrt/sources.py +++ b/pinballrt/sources.py @@ -51,6 +51,7 @@ def emit(self, nphotons, distance_unit, wavelength="random", simulation="thermal phi = 2*np.pi*np.random.rand(nphotons) direction = cost[:,np.newaxis]*r_hat + (sint*np.cos(phi))[:,np.newaxis]*phi_hat + (sint*np.sin(phi))[:,np.newaxis]*theta_hat + direction_frame = cost[:,np.newaxis]*r_hat + (sint*np.cos(phi))[:,np.newaxis]*phi_hat + (sint*np.sin(phi))[:,np.newaxis]*theta_hat if wavelength == "random": t1 = time.time() @@ -69,6 +70,7 @@ def emit(self, nphotons, distance_unit, wavelength="random", simulation="thermal photon_list = PhotonList() photon_list.position = wp.array(position, dtype=wp.vec3) photon_list.direction = wp.array(direction, dtype=wp.vec3) + photon_list.direction_frame = wp.array(direction_frame, dtype=wp.vec3) photon_list.frequency = wp.array(frequency, dtype=float) photon_list.energy = wp.array(photon_energy, dtype=float) @@ -91,6 +93,7 @@ def emit_rays(self, nu, distance_unit, ez, nrays, distance, device="cpu"): ray_list = PhotonList() ray_list.position = wp.array(position, dtype=wp.vec3) ray_list.direction = wp.array(direction, dtype=wp.vec3) + ray_list.direction_frame = wp.array(direction, dtype=wp.vec3) ray_list.indices = wp.zeros(position.shape, dtype=int) ray_list.intensity = wp.array2d(intensity, dtype=float) ray_list.tau_intensity = wp.array2d(tau_intensity, dtype=float) diff --git a/pinballrt/tests/test_E2E.py b/pinballrt/tests/test_E2E.py index 71ee42b..d730b54 100644 --- a/pinballrt/tests/test_E2E.py +++ b/pinballrt/tests/test_E2E.py @@ -5,10 +5,8 @@ from pinballrt.utils import calculate_Qvalue import astropy.units as u -import matplotlib.pyplot as plt import numpy as np import xarray as xr -import torch import os import pytest @@ -37,14 +35,21 @@ def test_E2E(grid_class, grid_kwargs, percentile, return_vals=False): # Set up the grid. model = Model(grid=grid_class, grid_kwargs=grid_kwargs) - density = np.ones(model.grid.shape)*1.0e-16 * u.g / u.cm**3 + density = np.ones(model.grid.shape)*1.0e-14 * u.g / u.cm**3 - model.add_density(density, d) + vx, vy, vz = np.meshgrid(0.5*(model.grid.grid.w1.numpy()[1:] + model.grid.grid.w1.numpy()[0:-1]), + 0.5*(model.grid.grid.w2.numpy()[1:] + model.grid.grid.w2.numpy()[0:-1]), + 0.5*(model.grid.grid.w3.numpy()[1:] + model.grid.grid.w3.numpy()[0:-1]), indexing='ij') + velocity = np.concatenate((vx[np.newaxis], vy[np.newaxis], vz[np.newaxis]), axis=0) * (-1.0 * u.km / u.s) + + model.set_physical_properties(density=density, dust="yso.dst", amax=1.0*u.mm, p=3.5, gases=['co.dat'], + abundances=[1.0e-4], microturbulence=0.2 * u.km / u.s, velocity=velocity) model.add_star(star) model.thermal_mc(nphotons=1000000, use_ml_step=False, Qthresh=1.02, Delthresh=1.02) - image = model.make_image(npix=256, pixel_size=0.2*u.arcsec, lam=np.array([1., 1000.])*u.micron, incl=45.*u.degree, pa=45.*u.degree, distance=1.*u.pc, nphotons=1000000) + image = model.make_image(npix=256, pixel_size=0.2*u.arcsec, channels=np.array([1., 1000.])*u.micron, + incl=45.*u.degree, pa=45.*u.degree, distance=1.*u.pc, nphotons=1000000, include_gas=False) # Do the checks. @@ -60,7 +65,7 @@ def test_E2E(grid_class, grid_kwargs, percentile, return_vals=False): base_image = xr.open_dataset(os.path.join(os.path.dirname(__file__), f"data/{grid_class.__name__}_E2E_image.nc")) Q = calculate_Qvalue(image.intensity, base_image.intensity, percentile=100.0, clip=0.1) - assert Q < 1.02, f"Image difference exceeds tolerance: {Q}" + assert Q < 1.03, f"Image difference exceeds tolerance: {Q}" else: return model.grid.grid.temperature.numpy(), model.grid.scattering.numpy(), image