{ "cells": [ { "cell_type": "markdown", "id": "a9441644", "metadata": {}, "source": [ "# Working with snapshots\n", "\n", "The core of `snapAnalysis` is the `snapshot` class, which handles data from one particle type of a single snapshot. In this notebook, we demonstrate a few common workflows for loading and working with simulation particle data. " ] }, { "cell_type": "code", "execution_count": 1, "id": "c69fe406", "metadata": {}, "outputs": [], "source": [ "from snapanalysis.snap import snapshot\n", "import astropy.units as u # for unit handling" ] }, { "cell_type": "markdown", "id": "3da6552f", "metadata": {}, "source": [ "## The snapshot object\n", "\n", "For this tutorial, we'll use one of the small snapshots included with `snapAnalysis`'s test suite, which is a low-resolution isolated dark matter halo. To instantiate a `snapshot` object, provide the filename and the particle type you want to work with:" ] }, { "cell_type": "code", "execution_count": 2, "id": "3538f772", "metadata": {}, "outputs": [], "source": [ "s = snapshot(\"../tests/example_snaps/cdm_snaps/snap_000.hdf5\", 1)" ] }, { "cell_type": "markdown", "id": "767e0ace", "metadata": {}, "source": [ "Instantiating a `snapshot` object loads the snapshot's metadata from its header, including time, units (handled via astropy), number of particles, etc. " ] }, { "cell_type": "code", "execution_count": 3, "id": "b9204247", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.0 Gyr 0.9777923542981724 Gyr\n", "43021.931322710916 9.99703e-11 km2 kpc / (solMass s2) 1.000000135629411 kpc 10002967845.365513 solMass 1.0 km / s\n", "10000\n", "0.0 kpc\n" ] } ], "source": [ "print(s.time, s.time_unit)\n", "print(s.G, s.length_unit, s.mass_unit, s.vel_unit)\n", "print(s.N)\n", "print(s.box_size)" ] }, { "cell_type": "markdown", "id": "d477d123", "metadata": {}, "source": [ "There are also convenience methods for printing the entire snapshot structure, or the entire header." ] }, { "cell_type": "code", "execution_count": 4, "id": "3db76b37", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "File structure of ../tests/example_snaps/cdm_snaps/snap_000.hdf5\n", "\n", "GROUP: Config\n", "KEYS:\n", " ALLOW_HDF5_COMPRESSION\n", " DOUBLEPRECISION\n", " DOUBLEPRECISION_FFTW\n", " EVALPOTENTIAL\n", " IMPOSE_PINNING\n", " MAX_NUMBER_OF_RANKS_WITH_SHARED_MEMORY\n", " NSOFTCLASSES\n", " NTYPES\n", " NUMBER_OF_MPI_LISTENERS_PER_NODE\n", " OUTPUT_ACCELERATION\n", " OUTPUT_IN_DOUBLEPRECISION\n", " OUTPUT_POTENTIAL\n", " OUTPUT_TIMESTEP\n", " POSITIONS_IN_64BIT\n", " RANDOMIZE_DOMAINCENTER\n", " SELFGRAVITY\n", "DATA SETS:\n", "\n", "GROUP: Header\n", "KEYS:\n", " BoxSize\n", " Git_commit\n", " Git_date\n", " MassTable\n", " NumFilesPerSnapshot\n", " NumPart_ThisFile\n", " NumPart_Total\n", " Redshift\n", " Time\n", "DATA SETS:\n", "\n", "GROUP: Parameters\n", "KEYS:\n", " ActivePartFracForNewDomainDecomp\n", " ArtBulkViscConst\n", " BoxSize\n", " ComovingIntegrationOn\n", " CourantFac\n", " CpuTimeBetRestartFile\n", " DesNumNgb\n", " ErrTolForceAcc\n", " ErrTolIntAccuracy\n", " ErrTolTheta\n", " ErrTolThetaMax\n", " GravityConstantInternal\n", " Hubble\n", " HubbleParam\n", " ICFormat\n", " InitCondFile\n", " InitGasTemp\n", " MaxFilesWithConcurrentIO\n", " MaxMemSize\n", " MaxNumNgbDeviation\n", " MaxSizeTimestep\n", " MinEgySpec\n", " MinSizeTimestep\n", " NumFilesPerSnapshot\n", " Omega0\n", " OmegaBaryon\n", " OmegaLambda\n", " OutputDir\n", " OutputListFilename\n", " OutputListOn\n", " SnapFormat\n", " SnapshotFileBase\n", " SofteningClassOfPartType0\n", " SofteningClassOfPartType1\n", " SofteningClassOfPartType2\n", " SofteningClassOfPartType3\n", " SofteningClassOfPartType4\n", " SofteningClassOfPartType5\n", " SofteningComovingClass0\n", " SofteningComovingClass1\n", " SofteningComovingClass2\n", " SofteningComovingClass3\n", " SofteningComovingClass4\n", " SofteningComovingClass5\n", " SofteningMaxPhysClass0\n", " SofteningMaxPhysClass1\n", " SofteningMaxPhysClass2\n", " SofteningMaxPhysClass3\n", " SofteningMaxPhysClass4\n", " SofteningMaxPhysClass5\n", " TimeBegin\n", " TimeBetSnapshot\n", " TimeBetStatistics\n", " TimeLimitCPU\n", " TimeMax\n", " TimeOfFirstSnapshot\n", " TopNodeFactor\n", " TypeOfOpeningCriterion\n", " UnitLength_in_cm\n", " UnitMass_in_g\n", " UnitVelocity_in_cm_per_s\n", "DATA SETS:\n", "\n", "GROUP: PartType1\n", "KEYS:\n", "DATA SETS:\n", " Acceleration\n", " Coordinates\n", " ParticleIDs\n", " Potential\n", " TimeStep\n", " Velocities\n", "\n", "GROUP: PartType5\n", "KEYS:\n", "DATA SETS:\n", " Acceleration\n", " Coordinates\n", " ParticleIDs\n", " Potential\n", " TimeStep\n", " Velocities\n" ] } ], "source": [ "s.print_structure()" ] }, { "cell_type": "code", "execution_count": 5, "id": "cc54cb0c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "___________________________________________\n", "Header of ../tests/example_snaps/cdm_snaps/snap_000.hdf5\n", "\n", "KEYS\n", "--------\n", "BoxSize: 0.0\n", "Git_commit: b'01e6b1567c93fe1cfaffd499aa55151db2ed4208'\n", "Git_date: b'Tue Mar 2 13:22:03 2021 +0100'\n", "MassTable: [0.00000000e+00 1.79999643e-03 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 1.79999645e-07]\n", "NumFilesPerSnapshot: 1\n", "NumPart_ThisFile: [ 0 10000 0 0 0 1]\n", "NumPart_Total: [ 0 10000 0 0 0 1]\n", "Redshift: 0.0\n", "Time: 0.0\n", "\n", "DATA SETS\n", "--------\n" ] } ], "source": [ "s.print_header()" ] }, { "cell_type": "markdown", "id": "dfd8642d", "metadata": {}, "source": [ "## Loading particle data\n", "\n", "`snapshot` objects do NOT automatically load particle data. Instead, the user can specify only which fields they want to work with, to save memory. To load (a) data field(s), use" ] }, { "cell_type": "code", "execution_count": 6, "id": "37f335d7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[-34.63337414 -66.0694746 172.63700156]\n", " [-36.79456065 -60.03053335 153.69383798]\n", " [-30.34222824 -17.90112547 161.11934558]\n", " ...\n", " [ 16.34056885 -11.39983141 171.08065064]\n", " [ 18.99446364 -65.16006591 128.64642553]\n", " [ 18.09111078 -53.39961967 171.65064279]] kpc\n", "[18005306.43454713 18005306.43454713 18005306.43454713 ...\n", " 18005306.43454713 18005306.43454713 18005306.43454713] solMass\n" ] } ], "source": [ "s.load_particle_data([\"Coordinates\", \"Velocities\"])\n", "\n", "# once loaded, the data is stored in a dictionary, indexed by field:\n", "print(s.data_fields[\"Coordinates\"])\n", "\n", "# you can load masses, positions, and velocities in a single line with\n", "s.read_all()\n", "print(s.data_fields[\"Masses\"])" ] }, { "cell_type": "markdown", "id": "0bcb4af7", "metadata": {}, "source": [ "Note that the particle loader will only load a field if it is not already loaded and no filters have been applied to the dataset." ] }, { "cell_type": "markdown", "id": "0c8a8293", "metadata": {}, "source": [ "## Centering and Rotations\n", "\n", "A common need when working with snapshots is to move to the center-of-mass (COM) frame of a galaxy, in which the z-axis is aligned with the galaxy's angular momentum vector. These can be accomplished quickly via" ] }, { "cell_type": "code", "execution_count": 7, "id": "693d5531", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/haydenfoote/Documents/gradschool/research/code/snapAnalysis/src/snapanalysis/snap.py:468: UserWarning: COM: Minimum number of particles reached before COM converged.\n", " warnings.warn(\n" ] } ], "source": [ "s.find_and_apply_center()\n", "s.align_angular_momentum()" ] }, { "cell_type": "markdown", "id": "ec9a0f38", "metadata": {}, "source": [ "You can also obtain the center of mass without applying the transformation with the following. Note these should be close to zero since we've already applied the centering in the previous cell." ] }, { "cell_type": "code", "execution_count": 8, "id": "cedcc6dc", "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$[-2.0729946 \\times 10^{-16},~3.9551695 \\times 10^{-16},~-8.500145 \\times 10^{-16}] \\; \\mathrm{kpc}$" ], "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "com = s.find_position_center()\n", "com" ] }, { "cell_type": "markdown", "id": "a77f3ec0", "metadata": {}, "source": [ "You'll notice the previous two cells gave a warning that the COM did not converge. The centering uses the shrinking-sphere alogrithm described in [Power et al. 2003](https://ui.adsabs.harvard.edu/abs/2003MNRAS.338...14P), in which the COM is computed using successively smaller spheres, stopping when the postiion converges. The test snapshots include too few particles to effectively apply this method using its default parameters, though you can adjust the shrinking rate of the sphere, initial size of the sphere, and other parameters using the kwargs:" ] }, { "cell_type": "code", "execution_count": 9, "id": "99312863", "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$[0.10042075,~0.22558081,~0.12055743] \\; \\mathrm{kpc}$" ], "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s.find_position_center(vol_dec=1.1, r_start=15.*u.kpc, guess=[0,0,0]*u.kpc)" ] }, { "cell_type": "code", "execution_count": 10, "id": "4133f806", "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$[-3.400058 \\times 10^{-15},~-5.2735594 \\times 10^{-16},~1.3908145 \\times 10^{-15}] \\; \\mathrm{\\frac{km}{s}}$" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# find velocity center requires the position center, as it finds the average velocity of particles within a small radius of the position center. \n", "com_v = s.find_velocity_center(com, r_max=15.*u.kpc)\n", "com_v" ] }, { "cell_type": "code", "execution_count": 11, "id": "619e8fce", "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$[5.7378153 \\times 10^{-15},~-1.2200395 \\times 10^{-15},~1] \\; \\mathrm{}$" ], "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# angular momentum should be close to [0,0,1]\n", "s.find_angular_momentum_direction()" ] }, { "cell_type": "markdown", "id": "ce1770b2", "metadata": {}, "source": [ "## Filtering particles\n", "\n", "The `select_particles` method allows one to keep a subset of particles according to various criteria. It supports ID ranges and boolean masks. For example, to keep all particles within 15 kpc of the center, use" ] }, { "cell_type": "code", "execution_count": 12, "id": "1972426b", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "r = np.sqrt(np.sum(s.data_fields[\"Coordinates\"]**2, axis=1)) # radial coordinate of all particles\n", "s.select_particles(r < 15*u.kpc)" ] }, { "cell_type": "markdown", "id": "4a1a2bbe", "metadata": {}, "source": [ "## Plotting\n", "`snapAnalysis` includes a visualization module with a few basic analysis and plotting tools. If you want to plot something that isn't included out-of-the-box, it should be straightforward to make your own plotting functions using the extracted particle data. If you make your own visualization functions, please consider [contributing them](/CONTRIBUTING.md) through a pull request, or [open an issue](https://github.com/hfoote/snapAnalysis/issues) so that we can implement the plot you want! \n", "\n", "Note that by default, the visualization functions will show a rudimentary version of the plot, but will also return the arrays used to make the plot so you can make them pretty in your own style. For example, let's plot the spherically-averaged density profile of the particles we've selected. " ] }, { "cell_type": "code", "execution_count": 13, "id": "472104a2", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAHACAYAAABTSTnVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGA0lEQVR4nO3deXxTdb7/8Xeatume0tKFlpZ9EyiyDyJuoILoiAsqoiLoXBdcGK/jMr+rjlcUxRkuV2cGtyvqVURnRhzlugEDiMoOZRMKCEiBtiylSde0TfL7ozS0toXSpD1J+3o+HudRcnJy8iHDNG+/q8ntdrsFAAAQYIKMLgAAAKApCDEAACAgEWIAAEBAIsQAAICARIgBAAABiRADAAACEiEGAAAEJEIMAAAISIQYAAAQkAgxAAAgILWJEPPtt9/qmmuuUUpKikwmkz799NNzev0f/vAHmUymOkdkZGTzFAwAAM6qTYSY4uJiDRgwQH/5y1+a9PpHH31UOTk5tY7zzjtPEydO9HGlAACgsdpEiBk3bpxmzpyp6667rt7nHQ6HHn30UaWmpioyMlLDhw/XihUrPM9HRUUpOTnZc+Tl5enHH3/UXXfd1UJ/AwAA8EttIsSczQMPPKDVq1dr4cKF2rp1qyZOnKixY8dqz5499V7/1ltvqWfPnho1alQLVwoAAKq1+RBz8OBBzZ8/X3/72980atQodevWTY8++qguvPBCzZ8/v871ZWVl+uCDD2iFAQDAYMFGF2C0bdu2yel0qmfPnrXOOxwOxcfH17l+0aJFKiws1JQpU1qqRAAAUI82H2KKiopkNpu1ceNGmc3mWs9FRUXVuf6tt97S1VdfraSkpJYqEQAA1KPNh5iBAwfK6XTq6NGjZx3jsn//fi1fvlyfffZZC1UHAAAa0iZCTFFRkfbu3et5vH//fmVmZiouLk49e/bU5MmTdccdd+hPf/qTBg4cqGPHjmnZsmXKyMjQ+PHjPa97++231aFDB40bN86IvwYAAKjB5Ha73UYX0dxWrFihSy+9tM75KVOm6J133lFFRYVmzpyp9957T4cPH1b79u31q1/9Ss8++6z69+8vSXK5XOrUqZPuuOMOPf/88y39VwAAAL/QJkIMAABofdr8FGsAABCYCDEAACAgteqBvS6XS0eOHFF0dLRMJpPR5QAAgEZwu90qLCxUSkqKgoIabm9p1SHmyJEjSktLM7oMAADQBNnZ2erYsWODz7fqEBMdHS2p6kOIiYkxuBoAANAYdrtdaWlpnu/xhrTqEFPdhRQTE0OIAQAgwJxtKAgDewEAQEAixAAAgIBEiAEAAAGJEAMAAAISIQYAAAQkQgwAAAhIhBgAABCQCDEAACAgEWIAAEBA8tsQ43Q69dRTT6lLly4KDw9Xt27d9Nxzz8ntdhtdGgAA8AN+u+3ASy+9pHnz5undd99V3759tWHDBk2dOlVWq1UPPfSQ0eUBAACD+W2I+eGHH3Tttddq/PjxkqTOnTvrww8/1Lp16wyuDAAA+AO/7U664IILtGzZMu3evVuStGXLFn333XcaN26cwZWdu7IKJ91gAAD4mN+2xDzxxBOy2+3q3bu3zGaznE6nnn/+eU2ePLnB1zgcDjkcDs9ju93eEqWe0ZGCUo2Zs1Jj+yZrzs3nG10OAACtht+2xHz88cf64IMPtGDBAm3atEnvvvuu/vjHP+rdd99t8DWzZs2S1Wr1HGlpaS1Ycf02/nxSJeVOrTuQb3QpAAC0Kia3n/ZzpKWl6YknntD06dM952bOnKn3339fu3btqvc19bXEpKWlyWazKSYmptlrrs9fV+zV7K+yFG0J1rZnrzSkBgAAAondbpfVaj3r97ffdieVlJQoKKh2Q5HZbJbL5WrwNRaLRRaLpblLOyeHTpZKkgodlap0uhRs9tvGLwAAAorfhphrrrlGzz//vNLT09W3b19t3rxZc+bM0bRp04wu7Zxk55d4/mwvq1RcZKiB1QAA0Hr4bYh59dVX9dRTT+n+++/X0aNHlZKSonvuuUdPP/200aWdk8OnWmIkyVZaQYgBAMBH/DbEREdHa+7cuZo7d67RpTSZy+X2dCdJUkFJuaRI4woCAKAVYYBGMzpa6FC58/QYnoLSCgOrAQCgdSHENKNDJ0tqPbYTYgAA8BlCTDPK/kWIKSghxAAA4CuEmGZ0KL+01mMbLTEAAPgMIaYZVbfEBJmqHtMSAwCA7xBimlH2qZaY7olRkmiJAQDAlwgxzehQQVVLTL9UqyTJVlpuZDkAALQqhJhmUul06UhBmSSpX0pViKE7CQAA3yHENJMcW5mcLrdCzUHqmRQtie4kAAB8iRDTTKpX6k1tF67YiBBJLHYHAIAvEWKaSfXMpI41QoyttEJut9vIsgAAaDUIMc2kuiWmY7sIWcOrQkx5pUtlFa4zvQwAADQSIaaZHMqvaolJiwtXlCVY5lOLxRQwQwkAAJ8gxDST091JETKZTIoNP92lBAAAvEeIaSbV3Ulp7cIlydOlxDRrAAB8gxDTDByVTuXaq9aISYuLkCRZI2iJAQDAlwgxzSCnoExutxQeYlZ8ZKik0y0xNlpiAADwCUJMM6g5vdpkqhrQy5gYAAB8ixDTDKo3fux4ajyMJMVGVLXIMDsJAADfIMQ0g0Mnq6dXR3jOxTCwFwAAnyLENINsz8yk0yGG7iQAAHyLENMMDtUYE1PNSogBAMCnCDHNoHpMTM3upFimWAMA4FOEGB8rLXfqeJFDUv0tMYyJAQDANwgxPna4oKorKdoS7AkuEi0xAAD4GiHGxzzTq+MiPGvESJI1vGqKtb2sQk6X25DaAABoTQgxPlbfoF7pdHeS2y0VltEaAwCAtwgxPlbf9GpJCg0OUkSoWRJdSgAA+AIhxsey86sXuguv8xyDewEA8B1CjI8dOlm95UBEnedYKwYAAN8hxPhY9slGtMQQYgAA8BohxocKyyo8XUX1tcR4plmXsAkkAADe8tsQ07lzZ5lMpjrH9OnTjS6tQdVdSe0iQhRlCa7zfOypadZ0JwEA4L2637R+Yv369XI6nZ7H27dv1+WXX66JEycaWNWZnR7UW7cVRpKsEQzsBQDAV/w2xCQkJNR6/OKLL6pbt266+OKLDaro7E4P6q07HkZiYC8AAL7ktyGmpvLycr3//vt65JFHaq2C+0sOh0MOh8Pz2G63t0R5Hp5BvfWMh5EY2AsAgC/57ZiYmj799FMVFBTozjvvPON1s2bNktVq9RxpaWktU+Apni0HGmiJYf8kAAB8JyBCzP/8z/9o3LhxSklJOeN1Tz75pGw2m+fIzs5uoQqreLYcaGhMTHV3EmNiAADwmt93J/38889aunSpPvnkk7Nea7FYZLFYWqCqutxut2dMTEPdSdWzkwpKmWINAIC3/L4lZv78+UpMTNT48eONLuWMbKUVKnJUSqI7CQCAluDXIcblcmn+/PmaMmWKgoP9u9GoejxMQrRFYSHmeq+JOdWdVFbhUlmFs95rAABA4/h1Mli6dKkOHjyoadOmGV1KLQ99uFm59jKFh5irjlCzThRXdRGlNdAKI0nRlmAFmSSXW7KXVjQYdgAAwNn5dYi54oor5Ha7jS6jjq2HCnTgREm9z3VNiGrwdUFBJsWEh6igpEIFpRVKjAlrrhIBAGj1/DrE+Kvnr+uvgpIKlVY4q47ySpWWu+SWWxOHnHlad+ypEMO4GAAAvEOIaYKR3ds3+bWeBe+YZg0AgFf8emBva2SNODXNmp2sAQDwCiGmhcWyfxIAAD5BiGlhbAIJAIBvEGJaGAveAQDgG4SYFsbAXgAAfIMQ08I8IYaWGAAAvEKIaWGxp2Yn0Z0EAIB3CDEtzDOwlynWAAB4hRDTwhjYCwCAbxBiWljNKdYul//tCwUAQKAgxLSw6hDjcktF5ZUGVwMAQOAixLSwsBCzLMFVH7uNadYAADQZIcYA1eNiWCsGAICmI8QYIDacadYAAHiLEGOA0wveMc0aAICmIsQYwMo0awAAvEaIMQD7JwEA4D1CjAFiT4UYOy0xAAA0GSHGALTEAADgPUKMATxTrBnYCwBAkxFiDGBlJ2sAALxGiDEA3UkAAHiPEGMABvYCAOA9QowBTi92R4gBAKCpCDEGqB7YW1LuVHmly+BqAAAITIQYA0SHhXj+zOBeAACahhBjAHOQSTFhwZIkG9OsAQBoEkKMQWKZZg0AgFcIMQZhmjUAAN4hxBgklp2sAQDwil+HmMOHD+u2225TfHy8wsPD1b9/f23YsMHosnwihpYYAAC8Emx0AQ05efKkRo4cqUsvvVRffvmlEhIStGfPHrVr187o0nwilrViAADwit+GmJdeeklpaWmaP3++51yXLl0MrMi3qruTWLUXAICm8dvupM8++0xDhgzRxIkTlZiYqIEDB+rNN98842scDofsdnutw1+dHtjLFGsAAJrCb0PMvn37NG/ePPXo0UNff/217rvvPj300EN69913G3zNrFmzZLVaPUdaWloLVnxuYsOZYg0AgDdMbrfbbXQR9QkNDdWQIUP0ww8/eM499NBDWr9+vVavXl3vaxwOhxwOh+ex3W5XWlqabDabYmJimr3mc/HV9lzd+/5GDUyP1aL7RxpdDgAAfsNut8tqtZ71+9tvW2I6dOig8847r9a5Pn366ODBgw2+xmKxKCYmptbhr5hiDQCAd/w2xIwcOVJZWVm1zu3evVudOnUyqCLfqh4TY2OKNQAATeK3Iea3v/2t1qxZoxdeeEF79+7VggUL9MYbb2j69OlGl+YT1S0xBaUV8tMePQAA/JrfhpihQ4dq0aJF+vDDD9WvXz8999xzmjt3riZPnmx0aT5RPbDX6XKruNxpcDUAAAQev10nRpKuvvpqXX311UaX0SzCQoIUag5SudOlgpJyRVn8+n8KAAD8jt+2xLR2JpNJVgb3AgDQZIQYAzG4FwCApiPEGKh6/yRaYgAAOHeEGANZz7IJZHZ+ia79y/ea8vY6ZjABAPALjCY1UPWYmIJ6upMyswt097vrdbyoam+lI7YypcaGt2h9AAD4M1piDNTQ/klf78jVLW+s9gQYScrK9d/NLAEAMAIhxkCegb2lp8PK29/t173vb1RZhUuX9ErQ6N6JkqSs3CJDagQAwF8RYgxUc/8kp8utP3y2Q/+5+Ee53dKtw9P11h1DNKhTO0m0xAAA8EuMiTFQdUtMjq1M976/UUt+zJMkPTGut+65qKtMJpN6JUVLknblFhpWJwAA/ogQY6Dqgb2bDxZIkkKDgzTnpgG6OiPFc02v5KoQ89OxIlU4XQox03gGAIBEd5KhqltiJKldRIgW3D28VoCRpNTYcEWGmlXhdGv/8eKWLhEAAL9FiDFQt4QoWcND1LV9pD65f6SGdI6rc01QkEk9T7XGZNGlBACAB91JBrKGh2j1k5fJEmyWOcjU4HW9k6O1+WCBsnILdc2AFiwQAAA/RogxWETo2f8nYHAvAAB10Z0UADzdSXlMswYAoBohJgD0To6RJGXnl6rIUWlwNQAA+AdCTACIiwxVQrRFkrQnjy4lAAAkQkzA6M0MJQAAaiHEBAgG9wIAUBshJkCwVgwAALURYgKEpzspr1But9vgagAAMB4hJkD0SIyWySTlF5freFG50eUAAGA4QkyACA81q3N8pCS6lAAAkAgxAeX04F4WvQMAgBATQBjcCwDAaYSYAFJzcC8AAG0dISaA9DoVYnbnFcrlYoYSAKBtI8QEkM7xkbIEB6mswqWD+SVGlwMAgKGCm/Kizz777Jxfc/nllys8PLwpb4dTzEEm9UiK0vbDdu3KLVTn9pFGlwQAgGGaFGImTJhwTtebTCbt2bNHXbt2bcrboYaeSdHaftiurNxCje2XbHQ5AAAYpsndSbm5uXK5XI06IiIifFlzm3Z6cC/TrAEAbVuTQsyUKVPOqWvotttuU0xMzDm9xx/+8AeZTKZaR+/evc+11FanV3LV58g0awBAW9ek7qT58+ef0/Xz5s1rytuob9++Wrp0qedxcHCTym1VqltiDpwoUVmFU2EhZoMrAgDAGH6dCoKDg5WczLiPmhKjLYqNCFFBSYX2Hi1Sv1Sr0SUBAGCIc+5OOnnypPLz8yVJx44d0yeffKIdO3b4vDBJ2rNnj1JSUtS1a1dNnjxZBw8ebJb3CSQmk0k9k1i5FwCAcwoxb731lgYPHqwhQ4Zo3rx5uu6667Rs2TLdcssteuutt3xa2PDhw/XOO+/oq6++0rx587R//36NGjVKhYUNf3E7HA7Z7fZaR2vEyr0AAJxjd9Irr7yiHTt2qLS0VOnp6dq/f78SEhJks9l08cUX6+677/ZZYePGjfP8OSMjQ8OHD1enTp308ccf66677qr3NbNmzdKzzz7rsxr8VS/2UAIA4NxaYoKDgxUeHq64uDh1795dCQkJkiSr1SqTydQsBVaLjY1Vz549tXfv3gavefLJJ2Wz2TxHdnZ2s9ZklN6EGAAAzi3EmM1mlZWVSZJWrlzpOV9UVOTbqupRVFSkn376SR06dGjwGovFopiYmFpHa1Q9JibXXiZbSYXB1QAAYIxzCjFLly6VxWKRVNX6Uq2kpERvvPGGTwt79NFHtXLlSh04cEA//PCDrrvuOpnNZk2aNMmn7xOIosNClBpbtU7PrtzWOe4HAICzOacQU1+30Q8//KCjR49q6NChPi3s0KFDmjRpknr16qWbbrpJ8fHxWrNmjacLq63rxeBeAEAb5/U6MdOnT9cDDzygfv361Tr/008/KTExUdHR0U2678KFC70trVXrlRytf+06yrgYAECb1eS9k6plZWXpkksuqXN+6dKldP00Iwb3AgDaOq9DTExMjE6ePFnn/KhRo7RmzRpvb48G1OxOcrvdBlcDAEDL8zrEjB07Vn/84x/r3jgoSOXl5d7eHg3o2j5KwUEmFZZV6oitzOhyAABocV6HmOeee04rV67UDTfcoG3btkmSysrK9NJLLykjI8PrAlG/0OAgdU2IlCRlMUMJANAGeR1i0tLStGbNGpWWlmrAgAEKDw9XdHS0Pv/8c7388su+qBEN6JVctQ5OVm7zr9MDAIC/8cku1p06ddIXX3yhgwcPavPmzQoNDdXw4cMVFxfni9ujAb2To/X5FlpiAABtk09CjCQdPnxYZrNZ1157ra9uibPodWrl3l3MUAIAtEFedyd9//336tKli9LT05Wenq6kpCQ9/vjjrXYHaX9SPUPpp2NFqnC6DK4GAICW5XWIueeee9SnTx+tX79eWVlZevnll7V06VINGjRIhw8f9kWNaEBqbLgiQ82qcLq1/3ix0eUAANCivA4xP/30k+bOnatBgwape/fuuuOOO7RhwwYNHDhQM2bM8EGJaEhQkEk9WfQOANBGeR1i+vTpo6NHj9Y6ZzKZ9J//+Z/66quvvL09zoKVewEAbZXXIebOO+/Ugw8+qOzs7FrnbTabYmJivL09zqIng3sBAG2U17OTqruMevTooeuvv17nn3++nE6n3n//fc2ePdvb2+MsTm8/wEBqAEDb4nWIycnJUWZmprZs2aLMzEy988472rNnj0wmk2bPnq0vv/xSGRkZysjI0NixY31RM2rofWrBu+z8UhU5KhVl8dmseQAA/JrX33gmk0lXXnmlrrzySs+5srIybdu2zRNuPvvsM73wwgsqKCjw9u3wC3GRoUqItuhYoUM7c+wa2pkFBgEAbYPJ7eUWyBdddJGWL18us9lc57nKykoFBxvXMmC322W1Wlv9+JzpCzbp/7bmaFiXOC38za8UFGQyuiQAAJqssd/fXg/sjY2N1UMPPVTn/IkTJzRmzBhvb49GeGJsb0WEmrVuf77eX/uz0eUAANAivA4x7733npYsWaK3337bc27nzp0aNmyYIiMjvb09GiEtLkJPjOstSXrxy13Kzi8xuCIAAJqfT1pi/vGPf+h3v/ud1q1bp6+//lojRozQhAkT9Pnnn/uiRjTCbcM7aViXOJWUO/XEJ1vlZS8hAAB+r0kDVqqnUlcf/fv315///GddddVVKisr06uvvqqpU6f6ulacQVCQSbNvyNDY//5W3+89oYXrszVpWLrRZQEA0Gya1BLTrVs3rVq1Snfffbc6d+6s+Ph4vfnmm3K73br11ls1aNAgVVRU+LpWnEXn9pF69IpekqTn/2+njhSUGlwRAADNx+vZSYcPH1ZmZmatY9++fQoODlbv3r21ZcsWX9V6ztrK7KSanC63bnztB20+WKBLeyXo7TuHymRithIAIHA09vvb6/nPqampSk1N1fjx4z3nioqKPGvEoGWZg0x6+cYMXfXKd1qedUyfbDqsGwZ3NLosAAB8rskDe59++mlt3Lix3ueioqJ04YUXavr06U0uDE3XPTFaM8b0kCQ9+/kOHbWXGVwRAAC+1+QQc+jQIY0bN04dO3bUfffdpy+//FLl5eW+rA1e+LdRXdU/1Sp7WaX+49PtzFYCALQ6TQ4xb7/9tnJzc/Xhhx8qOjpaM2bMUPv27XXDDTfovffeU35+vi/rxDkKNgfp5YkZCjGb9M2PeVq8NcfokgAA8Cmv1okJCgrSqFGjNHv2bGVlZWnt2rUaPny4Xn/9daWkpOiiiy7SH//4Rx0+fNhX9eIc9E6O0fRLu0uSnvlsh04UOQyuCAAA3/F6sbua+vTpo8cee0zff/+9Dh48qDvvvFOrVq3Shx9+6Mu3wTm4/5Lu6p0crfzicj3z2Q6jywEAwGe8nmLdkOLiYu3evVtpaWlq3759c7zFWbXFKdb12X7Ypmv/8r2cLrdeu22wxvZLNrokAAAa1GIbQNbnhRde0PXXX69PPvlE999/v6ZNm6bSUhZeM0q/VKvuvbirJOk/Pt2ughIGYAMAAp/PQ8z8+fOVn5+vr7/+Ws8995w+/vhjjRs3TjNmzPD1W+EcPHhZD3VPjNLxIoeeW7zT6HIAAPCaz0PMRx99pCeeeEKSNG3aNOXl5WnixIlav369r98K5yAsxKyXbsiQJH2aeVhFjkqDKwIAwDs+DzEhISGe9WJGjhypiIgISZLL5fLqvi+++KJMJhMtOl4Y3Kmd0uLC5XS5tf4AU+ABAIHN5yHm7rvv1mOPPSan06m77rpL0dHR+q//+i+NGzeuyfdcv369Xn/9dWVkZPiw0rZpRNd4SdKan04YXAkAAN7xeYi59tprddVVV2n06NGaMmWKrrjiCtlsNs2cObNJ9ysqKtLkyZP15ptvql27dj6utu0Z0e1UiNlHiAEABDavN4Csz6233qpbb71V+fn5ateunVe7KE+fPl3jx4/XmDFjmhyEcNqvTrXEbDtsk72sQjFhIQZXBABA0zQ5xMycOVODBg3S4MGDlZSUVO81cXFxTS5MkhYuXKhNmzY1elCww+GQw3F6VVq73e7V+7dGHazh6hwfoQMnSrR+f75G96n/fzsAAPxdk0PM008/7WlhSU5O9gSa6p+pqaleFZadna2HH35YS5YsUVhYWKNeM2vWLD377LNevW9bMKJbvA6cKNHqn04QYgAAAavJK/YOHz5cOTk5mjp1qtq3b69NmzZp48aN2rVrl5xOpxISEjRo0CB98cUXTSrs008/1XXXXSez2ew553Q6ZTKZFBQUJIfDUes5qf6WmLS0tDa/Yu8v/TPzsB5emKl+qTFa/OAoo8sBAKCWxq7Y2+SWmLVr1+qdd97R73//ew0dOlRz5sxRt27d5HA4lJmZqU2bNmnz5s1Nvb1Gjx6tbdu21To3depU9e7dW48//nidACNJFotFFoulye/ZVlTPUNpxxC5bSYWsEYyLAQAEHq9mJ915553avXu3evXqpUGDBunJJ5+U0+nU8OHDdd999+mNN95o8r2jo6PVr1+/WkdkZKTi4+PVr18/b8pu8xJjwtQ1IVJut7R2P7OUAACByesp1lFRUZo9e7Y2bNig7du3q3v37nrvvfd8URuaUXVrzGqmWgMAApRP1omprKyUw+HQpEmT1LFjR02dOlX5+b5fEXbFihWaO3euz+/bFp1eL4aVewEAganJY2JefPFFbdu2Tdu2bdOuXbsUFhamjIwMDRs2TPfcc4+sVqsv64SPDe9SFWJ25th1srhc7SJDDa4IAIBz0+QQ8/vf/16dO3fWlClTNGnSJPXs2dOXdaGZJURb1CMxSnuOFmnt/hMa26+D0SUBAHBOmtydNGrUKJ04cULPPvusBg8erJEjR+rBBx/U/PnztWXLFjmdTl/WiWZQ3aW0mn2UAAABqMktMStXrpQk7dmzRxs3btSmTZu0adMmffDBByooKJDFYlH//v21bt06nxUL3/pV13i9t/pnxsUAAAKS13sn9ejRQz169NAtt9ziObd//35t2LDBq3Vi0Pyq91HKyivUiSKH4qNYYwcAEDiaZQPILl26qEuXLpo4cWJz3B4+EhcZqt7J0dqVW6g1+/I1PoNxMQCAwNGkMTFbt26Vy+Vq9PU7duxQZWVlU94KzexXnvVijhtcCQAA56ZJIWbgwIE6caLxg0FHjBihgwcPNuWt0MyqQwzjYgAAgaZJ3Ulut1tPPfWUIiIiGnV9eXl5U94GLeBXXeNkMkl7jxbpaGGZEqMbt2M4AABGa1KIueiii5SVldXo60eMGKHw8PCmvBWaWWxEqPokx+jHHLvW7MvXrwekGF0SAACN0qQQs2LFCh+XASON6BZ/KsScIMQAAAKGT/ZOQmDzjIth0TsAQAAhxEDDusQpyCTtO16sPHuZ0eUAANAohBjIGh6ivilVG3ayBQEAIFAQYiCpapaSJK3ZR4gBAAQGQgwk1dgMkhADAAgQPgkxFRUVys7OVlZWlvLzWTQtEA3tHCdzkEk/nyjRkYJSo8sBAOCsmhxiCgsLNW/ePF188cWKiYlR586d1adPHyUkJKhTp076zW9+o/Xr1/uyVjSj6LAQ9UtlXAwAIHA0KcTMmTNHnTt31vz58zVmzBh9+umnyszM1O7du7V69Wo988wzqqys1BVXXKGxY8dqz549vq4bzYBxMQCAQNKkELN+/Xp9++23Wr16tfr3768LLrhA/fv3V/fu3TVs2DBNmzZN8+fPV25uriZMmKBVq1b5um40gxFdGRcDAAgcJrfb7fbmBuHh4dqxY4e6du3qq5p8xm63y2q1ymazKSYmxuhy/F6xo1IDnv1GlS63Vj12qdLiGrc3FgAAvtTY72+vB/YOHTpU+/fv9/Y28AORlmBldDw1LobWGACAn/M6xDz44IP6/e9/r+zsbF/UA4N5tiAgxAAA/JzXIebmm2/W+vXr1bdvX91222166623tHHjRpWXl/uiPrSw6vVi1vx0Ql72NAIA0KyatIt1Tfv379eWLVuUmZmpLVu2aNasWTpw4ICCg4PVq1cvbd261Rd1ooUM7tROIWaTjtjKdDC/RJ3iI40uCQCAenkdYjp16qROnTrp17/+tedcYWGhMjMzCTABKCI0WAM6xmrDzye1+qcThBgAgN9qlm0HoqOjNWrUKE2fPr05bo9m5ulSYlwMAMCPNSnEHDx48JyuP3z4cFPeBgapuV4M42IAAP6qSSFm6NChuueee864rYDNZtObb76pfv366R//+EeTC0TLG9SpnULNQcqzO7T/eLHR5QAAUK8mjYn58ccf9fzzz+vyyy9XWFiYBg8erJSUFIWFhenkyZP68ccftWPHDg0aNEizZ8/WVVdd5eu60YzCQsw6Pz1W6/bna/W+E+qaEGV0SQAA1NGklpj4+HjNmTNHOTk5+vOf/6wePXro+PHjnj2SJk+erI0bN2r16tUEmAA1wrNeDLuSAwD8k1ezk8LDw3XjjTdq2LBh+vbbb2WxWDRw4EB1797dV/XBICO6xeu/l+3R6lPrxZhMJqNLAgCgFq9nJ73yyivq2rWr7r//ft19993q1auXhg0b5vX06nnz5ikjI0MxMTGKiYnRiBEj9OWXX3pbLhrp/LRYWYKDdLzIoZ+OFRldDgAAdXgdYp577jk98cQTKigokM1mU1ZWli688EKNGDFC3333XZPv27FjR7344ovauHGjNmzYoMsuu0zXXnutduzY4W3JaISwELMGpbeTJK2mSwkA4Id8sov1tm3b6nQhzZw5U59//rnWrl3rVYE1xcXF6eWXX9Zdd93VqOvZxdo7ryzbozlLdmt8/w76y+RBRpcDAGgjWmwX64yMDK1evbrO+ZtuuslnK/Y6nU4tXLhQxcXFGjFiRIPXORwO2e32WgearnrRu9X7TuhkMXthAQD8i9ch5k9/+pP+/d//XR999FGthdHWrl2rHj16eHXvbdu2KSoqShaLRffee68WLVqk8847r8HrZ82aJavV6jnS0tK8ev+2bkDHWMWEBSu/uFwXzV6uuUt3q7CswuiyAACQ5IPuJEn64osvdO+996qsrEznn3++ysvLtX37dv3v//6vxo0b1+T7lpeX6+DBg7LZbPr73/+ut956SytXrmwwyDgcDjkcDs9ju92utLQ0upO8sP5Avp7+5w7tzKlq1YqNCNE9F3XTlAs6KSLU6623AACoo7HdST4JMVJVgFi2bJlWrFihzZs3a/PmzTp58qTi4uLUv39/ZWRkaO7cuV69x5gxY9StWze9/vrrjbqeMTG+4XK59cX2HM1Zslv7jlWt4Ns+yqLpl3bTpGHpCgsxG1whAKA1afEQU5/s7GxlZmZ6Qs2iRYu8ut9ll12m9PR0vfPOO426nhDjW5VOl/6ZeURzl+1Wdn6pJKmDNUwPXtZDE4d0VIi5WfYTBQC0MX4RYrzx5JNPaty4cUpPT1dhYaEWLFigl156SV9//bUuv/zyRt2DENM8Kpwu/W3DIb36rz3KsZVJktLjIjRjTA9de36qzEEsjAcAaLqADzF33XWXli1bppycHFmtVmVkZOjxxx9vdICRCDHNrazCqQVrD+qvK/bqeFHV7KXuiVF65PKeGts3WUGEGQBAEwR8iPEFQkzLKCmv1Ls//KzXVv4kW2nV7KXzOsToP67uowu6tTe4OgBAoGmxdWKAiNBg3XdJN616/FI9PLqHoizB+jHHrjvnr1eevczo8gAArRQhBj4TExai317eU6seu1T9U60qr3Tp7xsPGV0WAKCVIsTA59pFhur2EZ0kSX/bkK1W3GMJADAQIQbNYnz/DooMNevAiRKt3c8GkgAA3yPEoFlEWoJ1zYAUSdLH67MNrgYA0BoRYtBsbhpatXfVF9tzZGfPJQCAjxFi0GwGpsWqR2KUyipc+izziNHlAABaGUIMmo3JZNLNp1pjPt5AlxIAwLcIMWhW1w1MVYjZpK2HbJ6dsAEA8AVCDJpVfJRFY/okSZI+YoAvAMCHCDFodtUDfD/NPCxHpdPgagAArQUhBs3uoh4J6mANU0FJhb7ZkWd0OQCAVoIQg2ZnDjLpxsEdJTHAFwDgO4QYtIiJg6u6lL7be1yHTpYYXA0AoDUgxKBFpMdH6IJu8XK7pb9tYFNIAID3CDFoMdVrxvx94yE5XWwKCQDwDiEGLebKvsmKCQvW4YJSfb/3uNHlAAACHCEGLSYsxKwJA1MlSR8xwBcA4CVCDFrUTUOqupSW7MjTyeJyg6sBAAQyQgxaVL9Uq/qmxKjc6dKizYeNLgcAEMAIMWhxNTeFdLsZ4AsAaBpCDFrctQNSFRocpF25hdp6yGZ0OQCAAEWIQYuzRoRoXL9kSazgCwBoOkIMDHHzqQG+n2UeUWk5m0ICAM4dIQaG+FXXeKXFhavQUakvt+cYXQ4AIAARYmCIoCCTbjq1n9JH6+lSAgCcO0IMDHPjkI4KMklr9+frwPFio8sBAAQYQgwM08Earot6JkhigC8A4NwRYmCo6gG+f994SJVOl8HVAAACCSEGhhrdJ0nxkaE6WujQyt3HjC4HABBACDEwVGhwkK6r3hSSAb4AgHNAiIHhqrch+NeuozpW6DC4GgBAoPDbEDNr1iwNHTpU0dHRSkxM1IQJE5SVlWV0WWgGPZKiNTA9VpUutz7ZdMjocgAAAcJvQ8zKlSs1ffp0rVmzRkuWLFFFRYWuuOIKFRczFbc1qh7g+xGbQgIAGsnkDpBvjGPHjikxMVErV67URRdd1KjX2O12Wa1W2Ww2xcTENHOF8EaRo1LDnl+qknKn/n7vCA3pHGd0SQAAgzT2+9tvW2J+yWar2u04Lq7hLzeHwyG73V7rQGCIsgRrfP8OkhjgCwBonIAIMS6XSzNmzNDIkSPVr1+/Bq+bNWuWrFar50hLS2vBKuGt6gG+/7ctR0WOSoOrAQD4u4AIMdOnT9f27du1cOHCM1735JNPymazeY7sbP6LPpAM7tROXRMiVVLu1OItR4wuBwDg5/w+xDzwwANavHixli9fro4dO57xWovFopiYmFoHAofJZKo1wBcAgDPx2xDjdrv1wAMPaNGiRfrXv/6lLl26GF0SWsD1gzoqOMikzQcLtCev0OhyAAB+zG9DzPTp0/X+++9rwYIFio6OVm5urnJzc1VaWmp0aWhGCdEWXdY7URIDfAEAZ+a3IWbevHmy2Wy65JJL1KFDB8/x0UcfGV0amln1AN+P1mfrmx25BlcDAPBXwUYX0JAAWb4GzeDingkalB6rTQcL9G//u1E3Demop6/pqyiL3/5zBQAYwG9bYtB2BZuD9OG//Ur3XNxVJpP08YZDGvff32rd/nyjSwMA+BFCDPySJdisJ8f10Uf/NkId24UrO79UN7+xWrO+3ClHpdPo8gAAfoAQA782rEucvnx4lG4a0lFut/T6yn269s/fa1cuqzEDQFtHiIHfiw4L0ewbB+iN2wcrPjJUu3IL9etXv9cb3/4kp4uxUwDQVhFiEDCu6Jusr2ZcpDF9ElXudOmFL3bp1jfX6NDJEqNLAwAYgBCDgJIQbdGbdwzRi9f3V0SoWWv352vs3FX6+8ZDzGgDgDaGEIOAYzKZdMuwdH358CgN7tRORY5KPfq3Lbr3/Y06UeQwujwAQAshxCBgdYqP1Mf3jNBjY3spxGzS1zvydOXcVVq2M8/o0gAALYAQg4BmDjLp/ku669PpI9UzKUrHixy6690NevKTbSp2VBpdHgCgGRFi0Cr0TbHqswcu1N0XdpHJJH247qCuemWVNv7MAnkA0FoRYtBqhIWY9R9Xn6cP7h6uFGuYfj5RoomvrdYfv85SpdNldHkAAB8jxKDVuaBbe33124t0/aBUudzSn5fv1b3vb1RpOSv9AkBrQohBqxQTFqI5N52vVycNlCU4SEt3HtWkN9cwewkAWhFCDFq1awakaMFvhis2IkSZ2QW6Yd4P+vlEsdFlAQB8gBCDVm9wpzj9474L1LFduA6cKNH1f/1BW7ILjC4LAOAlQgzahG4JUfrk/gvULzVGJ4rLdcsba/SvXawnAwCBjBCDNiMxOkwL/22ELuqZoNIKp37z3kYtXHfQ6LIAAE1EiEGbEmUJ1v9MGaKJgzvK6XLriU+2ac6S3ey7BAABiBCDNifEHKTZN2boodE9JEmvLNujx/6+VRWsJQMAAYUQgzbJZDLpkct7atb1/WUOMulvGw/p7nc3sFUBAAQQQgzatEnD0vXmHYMVHmLWyt3HdPMbq3W0sMzosgAAjUCIQZt3We8kLfy3Xyk+MlTbD9t1/V9/0E/HiowuCwBwFoQYQNKAtFj9474L1Dk+QodOluqGeT+weSQA+DlCDHBK5/aR+sd9F2hAWqwKSip065tr9dX2XKPLAgA0gBAD1BAfZdGHvxmu0b0T5ah06b4PNuq91QeMLgsAUA9CDPALEaHBev32wbp1eLrcbunpf+7Qi1/uksvFWjIA4E8IMUA9gs1Ben5CPz16RU9J0msrf9IjH2eqvJK1ZADAXxBigAaYTCY9cFkP/XHiAAUHmfRp5hHdOX+d7GUVRpcGABAhBjirGwd31Nt3DlVkqFk//HRCN722Wrk21pIBAKMRYoBGuKhngj66Z4QSoi3alVuoq19dpWf+uV3f7j4mR6XT6PIAoE0yuVvxznd2u11Wq1U2m00xMTFGl4NWIDu/RFPfWa+9R08vhhcZataoHgka3SdRl/ZOVPsoi4EVAkDga+z3t1+HmG+//VYvv/yyNm7cqJycHC1atEgTJkxo9OsJMWgOZRVOrdpzXMt25mnZrqM6VujwPGcySeenxWpMnySN7pOoXknRMplMBlYLAIGnsd/fwS1Y0zkrLi7WgAEDNG3aNF1//fVGlwNIksJCzLr8vCRdfl6SXC63th+xaenOo/rXrjxtP2zX5oMF2nywQC9/naXU2HCN7pOo0X2S9KuucbIEm40uHwBaDb9uianJZDLREgO/l2sr07JdefrXzqP6bu9xOWpMyY4INWtUj/Ya3TtJl/ZOVEI03U4AUJ9W0RIDBJpka5gmD++kycM7qbTcqR9+Ou5ppcmzO/T1jjx9vSNPJpM0oGOsRveuaqXp04FuJwA4V60qxDgcDjkcp8cn2O12A6tBWxceatboPkka3SdJbnc/7Thi19KdeVq286i2HbYpM7tAmdkF+tOS3UqxhumyU91OI7rGKyyEbicAOJtW1Z30hz/8Qc8++2yd83Qnwd/k2cv0r11HtWznUX2395jKKk53O4WHmDWiW7wyOlqV0dGq/qmxdD0BaFNaxeykmhoTYupriUlLSyPEwK+VVVR1Oy3bWRVqcu11F9JLjglT/45WZaRa1a+jVf1TrUzlBtBqtckxMRaLRRYLv9gRWMJCzLqsd5Iu652kmRPc2nHErvUH8rXtsE3bDtm091iRcu1lyv2xTEt+zPO8LjU2XP1SY5TRMVb9UquCTVxkqIF/EwBoWX4dYoqKirR3717P4/379yszM1NxcXFKT083sDKgeZhMJvVLtapfqtVzrthRqR9z7Np6yKZthwq07bBN+44X63BBqQ4XlOrrHaeDTcd24eqfaj3VahOrfqkxio0g2ABonfy6O2nFihW69NJL65yfMmWK3nnnnbO+ninWaK2KHJXacdimbYdt2nrIpu2ngk190uMiagQbq/qmWmUND2nhigGg8VrdmJimIMSgLbGXVWj74apAUx1sDpwoqffazvER6pd6euDwoE6xLMQHwG8QYkSIAWwlFdp+xOYZX7P1cIGy80vrXBcdFqwrzkvW1QM66MLu7RViZm9YAMYhxIgQA9TnZHF5rWCz4eeTtfZ/io0I0ZWnAs2IrvEKJtAAaGGEGBFigMZwudxafyBf/7ctR19sy9HxonLPc/GRoRrbL1lXZ6RoWJc4mYNYVRhA8yPEiBADnCuny621+07o8605+mp7jk6WVHieS4i26Kp+ybp6QIoGp7dTEIEGQDMhxIgQA3ijwunS6p9OaPHWI/pqe67sZZWe5zpYw3RV/w66OqODzk+LZd8nAD5FiBEhBvCV8kqXvtt7TIu35OibH/NU5DgdaFJjw3V1RgddnZGifqkxBBoAXiPEiBADNIeyCqe+3X1Mi7fmaOnOPJWUOz3PdYqP0Pj+VYGGnbkBNBUhRoQYoLmVlju1POuo/m9rjpbtyqu1kWXXhEhdnZGiqzM6qGdStIFVAgg0hBgRYoCWVOyo1LJdR7V4yxGt2H1M5ZWnA03PpChPoOmaEGVglQACASFGhBjAKIVlFVq6M0+Lt+To2z3HVOE8/WvmvA4xGp/RQddkpCg9PsLAKgH4K0KMCDGAP7CVVOjrH3P1f1tz9P3e46p0nf6Vk9HRqkt7JapPh2j1So5RelwEa9EAIMRIhBjA35wsLtdXO3K1eOsRrf7phFy/+O0TFhKkHonR6pkUrd7J0ep16kiMtjBIGGhDCDEixAD+7FihQ1/vyNXmgwXanVeo3XmFctQYR1NTbERI7WCTFK2eydGKCWM3bqA1IsSIEAMEEqfLrYP5JcrKtWtXblWo2ZVbqAPHi+u02FRLjQ1Xz6Qo9UqOUa/kKPVKilG3xEh25AYCHCFGhBigNSircGrv0SJl1Qg2u/MKlWMrq/d6c5BJXdpHelpseiVXteCktYtgqwQgQBBiRIgBWjNbSYWy8gqrjly7ducWaVeuvdb2CDWFh5jVMylKPT3BJkY9k6OUEMV4G8DfEGJEiAHaGrfbrVx7mbJyC6uOvKqfe44W1Vq3pqa4yFD1TIpS7+QY9UqOVqf4CCVGW9Q+yiJreAgBBzAAIUaEGABVKp0u/ZxfcjrcnAo4B04U60y/AUPNQWofFaqEU6EmIdpS588JURa1j7YoMtRM4AF8hBAjQgyAMystPzXe5lSX1K7cQh0pKNWxQkeD3VINCQ8xnwo4oTUCTpjaR4cq4RcBKCyEgcfAmTT2+zu4BWsCAL8SHmpW/45W9e9orfNcWYVTJ4rLdazQoeOFDh0rcuhYYdVxvPrPp36WlDtVWuHUwfwSHcwvOev7RocF127R+eXPU8/FR4UqxBzUHH91oFUgxABAPcJCzEqNDVdqbPhZry12VHqCzfEaYacq5JTrWNGpIFToULnTpcKyShWWVWrfseKz3jsuMvRUl9XpFp3E6DAlW8PUwVr1MzE6TKHBhB20PYQYAPBSpCVYkZZgdYqPPON1brdb9rLKWiGnZitPzQB0orhcTpdb+cXlyi8uV1Zew/c1maT2UZaqUBNTHW7CPSGngzVMSTFhdGOh1SHEAEALMZlMsoaHyBoeou6JZ97N2+Vy62RJeb0BJ9fuUJ6tTLn2MuXaylTudHme2ypbg/dsFxFSO9zEVIeccE/YibTwtYDAwb9WAPBDQUEmxUdZFB9lUe/khq9zu6taa3JsVYEmx16mXFup53GurUxHbKUqq3DpZEmFTpZUaGeOvcH7RYcFn27JiandbVUddmLCgpmJBb9AiAGAAGYynQ47/VLrDlCWTnVjlVYqx3463FT9rB12Ch2Vp8brFGl3XlGD7xkRaj4dbmJqd1tVh512Eayxg+ZHiAGAVs5kMskaESJrRIh6Jzc8XbWwrEJ59qqAU2/YsZepoKRCJeVO7TtWfMaByaHBQfWO0UmMtijSEqyIULMiQoMVaTErPNSsyNBghYeY2RoC54QQAwCQJEWHhSg6LETdE6MbvKa03Klce5lybKU1Qs6pn/aqc8eLylVe6dLPJ0r084mzTzmvKTzErEhLVcCpCjrmWqGn+nHd64JrPY4MDVaEper14SEsRNhaEWIAAI0WHmpWl/aR6tK+4ZlYjkqnjtodp1p0aoed40UOFZc7VVpeqeJyp0oclSqpcHpWTi6tqFpzRyr3Wc0mkxQRYlZ4jaATGXq6Bag67ESGBp8ORZZTrUM1ron8RbCyBAcRjgxGiAEA+JQl2Ky0uAilxUU06nq3262yCpeKyytV4nCqpKJSxQ6nSsorVVJe9bPY4VRpubPqmnKnih2VtR5Xn6v+c/Vrq+4vFZc7VVzu1PGGh/qcsyCTPKEnPNQsS3CQQoODFGo+9TP49DmL51xQjevMdc5ZTr3eElL7+dPnarw+OEiWYLPMbbgLjhADADCUyWRS+KkgoDPPPD8nLpdbZZVOTyAqdjhV+ouA5GkNqg5LtR5XhSRPWHI4PaszS5LLLRU5KlXkOLctKnzNHGSqEZyCaoUpS42wUztg/SIMmYNkCTHXer5mcKr5ml/eKyHaYtjK0oQYAECrFBRkOtVSEizJ4rP7Ol1ulVacDjvVQae80iWH01X1s7LqZ9XhVLnTJUeFS+U1n691rur1v7yu+lpHjfuUV7rkcv+iHtfpcNXSlvz2IvVIangcVXMixAAAcA7MQSZFWYIVZeDCgJXOGkHJWTPsOE+Hp/oCkdMlR4Wz1rlaoau+ezkbvq680mXolhd+H2L+8pe/6OWXX1Zubq4GDBigV199VcOGDTO6LAAADBNsDlKwOUiRvmtgCkh+vWPYRx99pEceeUTPPPOMNm3apAEDBujKK6/U0aNHjS4NAAAYzK9DzJw5c/Sb3/xGU6dO1XnnnafXXntNERERevvtt40uDQAAGMxvQ0x5ebk2btyoMWPGeM4FBQVpzJgxWr16db2vcTgcstvttQ4AANA6+W2IOX78uJxOp5KSkmqdT0pKUm5ubr2vmTVrlqxWq+dIS0triVIBAIAB/DbENMWTTz4pm83mObKzs40uCQAANBO/nZ3Uvn17mc1m5eXl1Tqfl5en5OT696W3WCyyWNr4UG0AANoIv22JCQ0N1eDBg7Vs2TLPOZfLpWXLlmnEiBEGVgYAAPyB37bESNIjjzyiKVOmaMiQIRo2bJjmzp2r4uJiTZ061ejSAACAwfw6xNx88806duyYnn76aeXm5ur888/XV199VWewLwAAaHtMbrfbffbLApPdbpfVapXNZlNMTIzR5QAAgEZo7Pe3346JAQAAOBNCDAAACEiEGAAAEJAIMQAAICD59ewkb1WPWWYPJQAAAkf19/bZ5h616hBTWFgoSeyhBABAACosLJTVam3w+VY9xdrlcunIkSOKjo6WyWQyuhyfsdvtSktLU3Z2NlPHa+BzqYvPpH58LnXxmdSPz6WulvhM3G63CgsLlZKSoqCghke+tOqWmKCgIHXs2NHoMppNTEwM/6eqB59LXXwm9eNzqYvPpH58LnU192dyphaYagzsBQAAAYkQAwAAAhIhJgBZLBY988wzslgsRpfiV/hc6uIzqR+fS118JvXjc6nLnz6TVj2wFwAAtF60xAAAgIBEiAEAAAGJEAMAAAISISaAzJo1S0OHDlV0dLQSExM1YcIEZWVlGV2WX3nxxRdlMpk0Y8YMo0sx3OHDh3XbbbcpPj5e4eHh6t+/vzZs2GB0WYZxOp166qmn1KVLF4WHh6tbt2567rnnzrqseWvz7bff6pprrlFKSopMJpM+/fTTWs+73W49/fTT6tChg8LDwzVmzBjt2bPHmGJb0Jk+l4qKCj3++OPq37+/IiMjlZKSojvuuENHjhwxruAWcLZ/KzXde++9MplMmjt3bovVJxFiAsrKlSs1ffp0rVmzRkuWLFFFRYWuuOIKFRcXG12aX1i/fr1ef/11ZWRkGF2K4U6ePKmRI0cqJCREX375pX788Uf96U9/Urt27YwuzTAvvfSS5s2bpz//+c/auXOnXnrpJc2ePVuvvvqq0aW1qOLiYg0YMEB/+ctf6n1+9uzZeuWVV/Taa69p7dq1ioyM1JVXXqmysrIWrrRlnelzKSkp0aZNm/TUU09p06ZN+uSTT5SVlaVf//rXBlTacs72b6XaokWLtGbNGqWkpLRQZTW4EbCOHj3qluReuXKl0aUYrrCw0N2jRw/3kiVL3BdffLH74YcfNrokQz3++OPuCy+80Ogy/Mr48ePd06ZNq3Xu+uuvd0+ePNmgiownyb1o0SLPY5fL5U5OTna//PLLnnMFBQVui8Xi/vDDDw2o0Bi//Fzqs27dOrck988//9wyRRmsoc/k0KFD7tTUVPf27dvdnTp1cv/Xf/1Xi9ZFS0wAs9lskqS4uDiDKzHe9OnTNX78eI0ZM8boUvzCZ599piFDhmjixIlKTEzUwIED9eabbxpdlqEuuOACLVu2TLt375YkbdmyRd99953GjRtncGX+Y//+/crNza31/yOr1arhw4dr9erVBlbmf2w2m0wmk2JjY40uxTAul0u33367fve736lv376G1NCq905qzVwul2bMmKGRI0eqX79+RpdjqIULF2rTpk1av3690aX4jX379mnevHl65JFH9Pvf/17r16/XQw89pNDQUE2ZMsXo8gzxxBNPyG63q3fv3jKbzXI6nXr++ec1efJko0vzG7m5uZKkpKSkWueTkpI8z0EqKyvT448/rkmTJrXp/ZReeuklBQcH66GHHjKsBkJMgJo+fbq2b9+u7777zuhSDJWdna2HH35YS5YsUVhYmNHl+A2Xy6UhQ4bohRdekCQNHDhQ27dv12uvvdZmQ8zHH3+sDz74QAsWLFDfvn2VmZmpGTNmKCUlpc1+Jjh3FRUVuummm+R2uzVv3jyjyzHMxo0b9d///d/atGmTTCaTYXXQnRSAHnjgAS1evFjLly9v1bt0N8bGjRt19OhRDRo0SMHBwQoODtbKlSv1yiuvKDg4WE6n0+gSDdGhQwedd955tc716dNHBw8eNKgi4/3ud7/TE088oVtuuUX9+/fX7bffrt/+9reaNWuW0aX5jeTkZElSXl5erfN5eXme59qy6gDz888/a8mSJW26FWbVqlU6evSo0tPTPb97f/75Z/37v/+7Onfu3GJ10BITQNxutx588EEtWrRIK1asUJcuXYwuyXCjR4/Wtm3bap2bOnWqevfurccff1xms9mgyow1cuTIOtPvd+/erU6dOhlUkfFKSkoUFFT7v9vMZrNcLpdBFfmfLl26KDk5WcuWLdP5558vSbLb7Vq7dq3uu+8+Y4szWHWA2bNnj5YvX674+HijSzLU7bffXmcM4pVXXqnbb79dU6dObbE6CDEBZPr06VqwYIH++c9/Kjo62tNHbbVaFR4ebnB1xoiOjq4zJigyMlLx8fFteqzQb3/7W11wwQV64YUXdNNNN2ndunV644039MYbbxhdmmGuueYaPf/880pPT1ffvn21efNmzZkzR9OmTTO6tBZVVFSkvXv3eh7v379fmZmZiouLU3p6umbMmKGZM2eqR48e6tKli5566imlpKRowoQJxhXdAs70uXTo0EE33nijNm3apMWLF8vpdHp+/8bFxSk0NNSospvV2f6t/DLIhYSEKDk5Wb169Wq5Ilt0LhS8IqneY/78+UaX5leYYl3l888/d/fr189tsVjcvXv3dr/xxhtGl2Qou93ufvjhh93p6enusLAwd9euXd3/7//9P7fD4TC6tBa1fPnyen+PTJkyxe12V02zfuqpp9xJSUlui8XiHj16tDsrK8vYolvAmT6X/fv3N/j7d/ny5UaX3mzO9m/ll4yYYs0u1gAAICAxsBcAAAQkQgwAAAhIhBgAABCQCDEAACAgEWIAAEBAIsQAAICARIgBAAABiRADAAACEiEGgN+48847ZTKZZDKZ9Omnn0qSDhw4IJPJpMzMzGZ733feecfzvjNmzGi29wHgW4QYAH5l7NixysnJ0bhx41rsPW+++Wbl5ORoxIgRLfaeALzHBpAAWlx5eXmDm+ZZLBYlJye3aD3h4eEKDw9vtRv5Aa0VLTEAmt0ll1yiBx54QDNmzFD79u115ZVXNvleTqdT06ZNU+/evXXw4EFJkslk0rx58zRu3DiFh4era9eu+vvf/17rdYcOHdKkSZMUFxenyMhIDRkyRGvXrvXq7wXAWIQYAC3i3XffVWhoqL7//nu99tprTbqHw+HQxIkTlZmZqVWrVik9Pd3z3FNPPaUbbrhBW7Zs0eTJk3XLLbdo586dkqSioiJdfPHFOnz4sD777DNt2bJFjz32mFwul0/+bgCMQXcSgBbRo0cPzZ49u8mvLyoq0vjx4+VwOLR8+XJZrdZaz0+cOFF33323JOm5557TkiVL9Oqrr+qvf/2rFixYoGPHjmn9+vWKi4uTJHXv3r3pfxkAfoGWGAAtYvDgwV69ftKkSSouLtY333xTJ8BIqjMod8SIEZ6WmMzMTA0cONATYAC0DoQYAC0iMjLSq9dfddVV2rp1q1avXn3Orw0PD/fqvQH4J0IMgIBw33336cUXX9Svf/1rrVy5ss7za9asqfO4T58+kqSMjAxlZmYqPz+/RWoF0DIIMQACxoMPPqiZM2fq6quv1nfffVfrub/97W96++23tXv3bj3zzDNat26dHnjgAUlVXVHJycmaMGGCvv/+e+3bt0//+Mc/mtSqA8B/EGIABJQZM2bo2Wef1VVXXaUffvjBc/7ZZ5/VwoULlZGRoffee08ffvihzjvvPElSaGiovvnmGyUmJuqqq65S//799eKLL8psNhv11wDgA8xOAtDsVqxY0eTXdu7cWW63u9a5Rx55RI888kitcykpKfrmm28avE+nTp3qrB0DILDREgPAryxevFhRUVFavHhxi73nBx98oKioKK1atarF3hOA90zuX/4nDgAY5OjRo7Lb7ZKkDh06NHpGk8lk0qJFizRhwoQmvW9hYaHy8vIkSbGxsWrfvn2T7gOgZRFiAABAQKI7CQAABCRCDAAACEiEGAAAEJAIMQAAICARYgAAQEAixAAAgIBEiAEAAAGJEAMAAAISIQYAAASk/w/Gp3mBOJhEcAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from snapanalysis import vis\n", "# this will be noisy since we're working with a small number of particles\n", "bins, density = vis.density_profile(s, rmin=1., rmax=15., nbins=20, log_bins=True)" ] }, { "cell_type": "code", "execution_count": 14, "id": "e00a54e9", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[None,\n", " None,\n", " Text(0.5, 0, 'Radius [kpc]'),\n", " Text(0, 0.5, 'Density [$M_\\\\odot / kpc^3$]')]" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAG1CAYAAAD+7yA/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+7klEQVR4nO3dd3hVVaKG8fckIXSilAFpgwULKiAYEAsDihQVRRisdMsMo4yI2BEbYkFRHCKjjqDYFRVRRxFQQUXBoSgoMuLYhqrSEQSS3D/WJciACDkn2UnO+3uePDkl2fuLc73nc+2114rl5ubmIkmSlGRSog4gSZIUBUuQJElKSpYgSZKUlCxBkiQpKVmCJElSUrIESZKkpGQJkiRJSckSJEmSklJa1AGKspycHJYsWULFihWJxWJRx5EkSXsgNzeXdevWUbNmTVJSfn28xxK0G0uWLKFOnTpRx5AkSfnw3XffUbt27V993xK0GxUrVgTCP8RKlSpFnEaSJO2JtWvXUqdOnbzP8V9jCdqNbZfAKlWqZAmSJKmY+a2pLE6MliRJSckSJEmSkpIlSJIkJSVLkCRJSkqWoF3IysqiQYMGZGZmRh1FkiQVkFhubm5u1CGKqrVr15KRkcGaNWu8O0ySpGJiTz+/HQmSJElJyRIkSZKSkiVIkiQlJUuQJElKSpYgSZKUlCxBkiQpKVmCtHe+/Rb694c77gBXV5AkFWPuIr8LWVlZZGVlkZ2dHXWUomPFChg6FEaNgs2bw2tlyoRCJElSMeRiibvhYonAmjVw991w772wYUN47fDD4dNPISUFJk6ENm2izShJ0i+4WKLis3EjDBsGBxwAQ4aEApSZCZMmwbx50KsX5OTA2WfDf/4TdVpJkvaaJUg72rIF/v53OOgguOoqWLkSDjsMXngBZswIoz6xWLgs1rx5eP+MM2D9+qiTS5K0VyxBCnJy4KmnQuHp2xeWLIHf/x4efTSM/HTuHMrPNmXKwIsvwn77wfz50LNnOIYkScWEJSjZ5ebCq6/CUUfB+efDl1/C734H998PCxeGcpOauuvfrVkzFKH09PD9ttsKN7skSXGwBCWzqVPh+OOhY0f45BPIyAjzf778Evr1g9Klf/sYxxwTLo0BDB4ML79csJklSUoQS1Aymj0b2reHVq1g+nQoWxauvjpMcL7+eqhQYe+O16dPKE0A3brBZ58lPLIkSYlmCUomn38OXbtC06bh1va0tDD/Z9GisPhh5cr5P/Y990Dr1mGC9BlnwKpVicstSVIBsAQlg2+/hQsuCOv7jBsXJjh36xbm/DzwQJjbE69SpeC558Jk6kWL4NxzwcUmJUlFmCWoJFuxAi6/HOrXh9Gjw91bp58OH38Mjz8e1gBKpKpVYfx4KFcujDRde21ijy9JUgK5bcYuFPi2GQsWwNdfb7/l/JffE/X4jTfCKs/b1u9p1Spse9GiRcH8Tds0bgxjxoRFFIcNC8/PO69gzylJUj64bcZuFNi2GVdcAcOHJ+54u9O0Kdx++/ZFDgvLddeF85YpA++9F3JIklQI9vTz25GgKNSuHdblge07sefm7v3j3b1frRoMHLjzIoeF5dZbw233r70GZ54J//pXWH9IkqQiwpGg3XAD1TitWRO21li4EE44ASZPDgsrSpJUgNxAVdHLyAiLJ1aqBO++C/37R51IkqQ8liAVrEMOCXuSbdt09aGHok4kSRJgCVJhOPXU7fuKXXopvP9+tHkkScISpMJyzTVw1lmwZQt06QLffRd1IklSkrMEqXDEYmHBxkaNYPnycMfYxo1Rp5IkJTFLkApP+fJhRekqVWDWLLj44u239EuSVMgsQSpc9eqF/ctSU+GJJ8Kq1pIkRcASpMLXqtX28nPllfDmm5HGkSQlJ0uQonHppdC7d9jU9Zxz4Msvo04kSUoylqBdyMrKokGDBmRmZkYdpeTatm7QMcfAqlVwxhmwbl3UqSRJScRtM3bDbTMKwZIlcPTRsHQpdOoEL7wAKXZzSVL+uYGqioeaNeGll6Bly3Dn2K23wo03xnfMrVthxYpQsJYuDd9/+Tg3F66/Ho49NiF/giSpeHIkaDccCSpEY8ZAnz7h8UsvhVGh//XLcvPLUvO/RWfFijDXaHfS0+Hhh6FHj4T/KZKkaDkSpOKld2+YOxfuvx+6d4fLLtu58CxfvufrCqWmQvXqYaRpv/3C922PX389FK2ePWHBgrClh5fgJCnpOBK0G44EFbItW6BdO3j77V//mV+Wm18WnP8tOtWqhZ/dlZwcGDx4+35mnTrB449DhQoJ/5MkSYVvTz+/LUG7YQmKwI8/hvk6sOuis7tys7eefBIuuAB+/jls5zFhAtStm5hjS5IiYwlKAEtQEvjggzAStGJFGGEaPz7cti9JKrb29PPbiRBKbi1awMyZ0LBhmHPUqhU8/XTUqSRJhcASJP3+9/Dee9CxY7g0dt55Yc7Qb91hJkkq1ixBEkDFiuGOsSuvDM9vvRXOPht++inaXJKkAmMJkrZJTYW77oLRo6FUqbDbfcuWsHhx1MkkSQXAEiT9r969YcoUqFIFZs2CzEz417+iTiVJSjBLkLQrJ5wQJkw3aBAWamzZEp5/PupUkqQEsgRJv+aAA8It9B06wMaNcNZZYa6Qq0pIUolgCdqFrKwsGjRoQGZmZtRRFLVKleCVV6B///B88OBw99jGjZHGkiTFz8USd8PFErWDhx6CSy4JG7lmZsLLL4dVrCVJRYqLJUqJdvHF8OabsO++8NFH0KwZzJkTdSpJUj5ZgqS90bo1zJgBhxwC//0vHH98WF9IklTsWIKkvVW/Pnz4IZx8clhMsXNnuP12J0xLUjFjCZLyY5994J//DHOEAK67Dnr0gE2bIo0lSdpzliApv9LSYOTI8JWaCk88ASeeGDZilSQVeZYgKV6XXAKvvw4ZGWFdoYYN4dJLYeLEsCGrJKlIsgRJiXDyyWGeUP36sGIFZGVB+/ZQtSp06QKPPgrffx91SknSL7hO0G64TpD22qZNMGlSWGDx1VfDlhvbxGLQogV07Bi+GjQIr0mSEmpPP78tQbthCVJccnJg9myYMCGUorlzd3x///23F6KWLSE9PZKYklTSWIISwBKkhPruuzA69Mor8NZbO84XqlQpXD7r2DHsVValSnQ5JamYswQlgCVIBWb9epg8eftlsxUrtr+XkgLHHRcK0emnh4UZJUl7zBKUAJYgFYqcHJg5MxSiV16BefN2fL9+/e2XzY4/PtyaL0n6VZagBLAEKRJff739stnbb8OWLdvf22cf6NcPbr7ZSdWS9CssQQlgCVLk1q0Lm7ZOmACvvQY//hhev/56GDIk2mySVES5i7xUElSsGNYZeuyxsBL13/4WXr/tNrj//mizSVIxZwmSiovU1LAS9bYRoMsug6eeijaTJBVjlqBdyMrKokGDBmRmZkYdRdrZddeFeUEAPXuG7TkkSXvNOUG74ZwgFVk5OdCtGzz9NJQvD1OmQPPmUaeSpCLBOUFSSZaSEvYja9sWNmyAU0+Fzz+POpUkFSuWIKm4Sk+HF16AZs3CXWNt28J//xt1KkkqNixBUnFWoUK4df6QQ8K2HG3bbr+NXpK0W5YgqbirWjWsJVSrFixYAKedFi6RSZJ2yxIklQR164a7xPbdFz78ELp23XGlaUnSTixBUklx+OHh0ljZsvD669CnT7iLTJK0S5YgqSRp0QLGjQsLKz7xBAwcCK6CIUm7ZAmSSppTToExY8Lje++Fu+6KNo8kFVGWIKkk6t4d7rknPL7mGhg9Oto8klQEWYKkkmrAALj66vD4oovCTvSSpDyWIKkku/126N07TJA++2yYNi3qRJJUZFiCpJIsFoOHHoLTT4dNm8L3Tz6JOpUkFQmWIKmkS0uDZ56B44+HNWugXTv46quoU0lS5CxBUjIoWzbMCTrySFi2LGyvsXx51KkkKVKWIClZ7LsvvPEG1KsHixZBhw6wdm3UqSQpMpYgKZnUrBn2GatWDebMgU6dwlwhSUpCliAp2dSvH7bVqFAB3n4bunWD7OyoU0lSobMEScmoaVMYPx7S0+GFF+CSS9xeQ1LSsQRJyeqkk8L+YrEYPPgg3HRT1IkkqVBZgqRk1rUrZGWFx7fcAiNHRptHkgpRWtQBJEWsb19YsSKMBP31rzBjBmRmhktmjRtD+fJRJ5SkAhHLzXUiwP/KysoiKyuL7Oxs/v3vf7NmzRoqVaoUdSyp4OTmwqWXwgMP7Ph6SgocemgoRNu+GjcOk6olqYhau3YtGRkZv/n5bQnajT39hyiVCLm5MGkSTJ8Os2aFr6VLd/65WGzHYtSkCRx1FFSsWPiZJWkXLEEJYAlS0lu6dHsh2va1ZMnOPxeLwcEH7zhidNRR4L83kiJgCUoAS5C0C8uWwezZOxaj//531z/7y2LUpEn47r9LkgqYJSgBLEHSHlqxYucRo+++2/nnypeHl16Ck08u/IySkoYlKAEsQVIcvv9+xxGjmTPDiFGNGjBvHlStGnVCSSWUJSgBLEFSAm3cCEcfDZ99Bp07w7hxYS6RJCXYnn5+u1iipMJRtmxYobpUKXjxRRg7NupEkpKcJUhS4TnqKLj55vC4Xz/4+utI40hKbpYgSYXrqqvguONg3Tro0cMd7CVFxhIkqXClpoZLYRUqwLvvwj33RJ1IUpKyBEkqfAccACNGhMeDBsHcuZHGkZScLEGSotG7N3TqBFu2QLdusGlT1IkkJRlLkKRoxGLw0EPwu9/Bp5/C9ddHnUhSkrEESYpOtWowenR4PHw4vPVWtHkkJRVLkKRonXoq/OlP4XHPnrB6daRxJCUPS5Ck6N19Nxx0UNhW45JLok4jKUlYgiRFr0KFsJp0aio89RQ880zUiSQlAUuQpKKhefPtk6P79oXFi6PNI6nEswRJKjoGDYLMzDAvqFcvyMmJOpGkEswSJKnoKFUKHn88bLY6eTKMHBl1IkklmCVIUtFyyCFhojTA1VfDZ59Fm0dSiWUJklT09O0L7duHVaS7dYPNm6NOJKkEsgRJKnpisbCIYuXKMGcO3Hxz1IkklUCWIElF0377hW01AO64A95/P9o8kkocS5CkoqtLl7CKdE4OdO8O69ZFnUhSCZKWn1+aMGHCXv/OySefTNmyZfNzOknJbMQIeOcd+Oor6N8fHnkk6kSSSohYbm5u7t7+UkrK3g0gxWIxvvjiCw444IC9PVWk1q5dS0ZGBmvWrKFSpUpRx5GS17Rp0KoV5ObCSy9Bp05RJ5JUhO3p53e+L4ctW7aMnJycPfoqV65cfk8jSdCyJVx5ZXh80UWwfHm0eSSVCPkqQT179tyrS1vdunVzJEVSfG65BRo1gh9+gAsuCKNCkhSHfF0OSxZeDpOKmPnzoWnTsG7Qgw/CxRdHnUhSEVTgl8MkqdAdcQTcfnt4fPnlsGhRtHkkFWt7XYJWrVrFypUrAfj+++958cUX+fTTTxMeTJJ2qX9/aN0afvoprCa9dWvUiSQVU3tVgv7xj3/QtGlTjj76aEaNGsWZZ57JlClTOOecc/jHP/5RUBkLXVZWFg0aNCAzMzPqKJL+V0oKPPooZGTAjBnbR4YkaS/t1Zyghg0bMmPGDDZu3EjdunX56quvqFatGmvWrOEPf/gDc+fOLcCohc85QVIR9uSTYSQoNRU++AD8jxZJ/69A5gSlpaVRtmxZKleuzEEHHUS1atUAyMjIIBaLxZdYkvbGeefBWWdBdnYoQz/9FHUiScXMXpWg1NRUNm3aBMDUqVPzXl+/fn1iU0nSb4nFYNQoqFkT/v1vuOqqqBNJKmb2qgRNnjyZ0qVLA2H0Z5uffvqJh7ZtdChJhaVy5TA/CCArC954I9I4koqXvSpBu7rsNX36dFasWOEkYknROPlk+Otfw+Pu3WHcOBdSlLRH4l4n6JJLLmHGjBk7vf7ll1+yzh2fJRWGO+7Yvpp0167Qti0sWBB1KklFXNwlaOHChbRq1Wqn1ydPnsy5554b7+El6beVLQvTp8PgwVC6NEyeDA0bhnlC/seYpF8RdwmqVKkSq1at2un1E044gQ8//DDew0vSnilXDm6+GT77DDp2DIsoDhsGhx4KTz/tJTJJO4m7BLVv356777575wOnpLB58+Z4Dy9Je+eAA2DCBHj1VTjwQFiyJNxO37p12HtMkv5f3CXo1ltvZerUqXTp0oV58+YBsGnTJu68804aNmwYd0BJypdTTw2l59Zbw+WyqVOhceOw7caaNVGnk1QExF2C6tSpw4cffsjGjRtp1KgRZcuWpWLFirzyyisMGzYsERklKX/KlIFBg8Ik6S5dwsKKI0bAwQfDY49BTk7UCSVFaK+2zfgt3377LXPmzCE9PZ3mzZtTuXLlRB06Em6bIZUwb74ZbqdfuDA8P/bYsL5Q48aRxpKUWHv6+Z2wErR48WIAatWqlYjDFQmWIKkE2rwZ7rsPbrkFNmwIG7L++c8wZAjsu2/U6SQlQIHsHbYr77//Pvvvvz9169albt26VK9enauvvpq1a9fGe2hJSrz09HDr/Oefw9lnh0tiDzwQLpH94x9eIpOSSNwl6E9/+hOHHXYYH330EQsXLmTYsGFMnjyZJk2a5I0OSVKRU7s2PPMMvPUWNGgQFlq86CI45hj46KOo00kqBHFfDitbtiwff/wxBx98cN5rubm5nHXWWQA8//zz8SWMkJfDpCSxZQuMHAk33hgWV4zF4MILYehQqFo16nSS9lKhXQ477LDDWLFixQ6vxWIxbrnlFt5wM0NJxUGpUnD55WE3+u7dw8KKDz8cLpGNGhXuKpNU4sRdgnr16kW/fv347rvvdnjd0RNJxU6NGjB2LLz7bth2Y9Uq+MtfIDMTPvgg6nSSEizuy2EpKaFHpaen07lzZxo3bkx2djZPPPEE1113Heeff35CgkbBy2FSEtu6FR58MKwztHp1eG3AALj77nC5TFKRVWi3yC9fvpy5c+fy8ccfM3fuXObOncsXX3xBLBbjsMMO48gjj6Rhw4Y0bNiQ9u3bx3OqQmcJksSKFXDddfDII+H5TTeFuUOSiqxCK0ErVqzgd7/73Q6vbdq0iXnz5u1QjubPn8/qbf81VUxYgiTl+fvfoW/f8PjBB+Hii6PNI+lXFVoJatmyJW+//Tapqak7vbd161bS0tLiOXykLEGSdjB4cNiLLCUFXnoJTj896kSSdqHQ7g7bZ599+Otf/7rT6z/++CNt2rSJ9/CSVHTcfDNccEFYUPHss2H69KgTSYpD3CVo7NixTJo0idGjR+e9tmDBApo1a0b58uXjPbwkFR2xWLgsdtppsGlT+L5gQdSpJOVTQkaCXnjhBa688kpmzpzJxIkTadGiBZ06deKVV15JREZJKjrS0sJK082bh1vo27cHV8eXiqV8TdjZdiv8tq8jjzySkSNHcsopp7Bp0yb+9re/0bt370RnlaSioXx5ePVVOO64sMBihw4wbRrss0/UySTthXyNBB144IG8++67XHjhhdSrV48qVarw8MMPk5uby3nnnUeTJk3YsmVLorNKUtFRtSpMnBgWWJw3Dzp1CpfIJBUbcd8dtnjx4rz1gbZ9/ec//yEtLY1DDz2Ujz/+OFFZC513h0n6TXPnQsuWYc+xP/4xXCrbxd2ykgpPod0ivyvr16/PWyPokksuSfThC40lSNIeeeutMDdoyxbo1w9GjHBVaSlCBX6L/ODBg5k1a9Yu36tQoQLHH398sS5AkrTHTjwRHn88PP7b3+Cuu6LNI2mP5LsE/fe//6VDhw7Url2bvn378vrrr7N58+ZEZpOk4uPss+Hee8Pja66Bxx6LNo+k35TvEjR69GiWLVvG008/TcWKFenfvz9Vq1alS5cujB07lpUrVyYypyQVff37w5VXhscXXABvvBFpHEm7l9A5QQsWLOCVV17h5ZdfZtasWTRr1ozTTz+dc889l1q1aiXqNIXGOUGS9lpODvToAU8+GW6lf/ttyMyMOpWUVCKdGA3w/fffM2HCBCZMmMAJJ5zAwIEDC+I0BcoSJClfNm8Oq0lPmgTVqsH770P9+lGnkpJG5CWoJLAEScq3deugVSuYPRv23x8++ACqV486lZQUCm0DVUnSLlSsCP/8JxxwAHz1FZxySihGkoqMAitBGzZsYM6cOfzwww8FdQpJKtqqVw+To6tVCyNCXbqES2WSioQCKUFDhw6lc+fOvPjii/zlL3+hT58+bNy4sSBOJUlFW/368NprUK5cmCPUp0+YPC0pcvnaQHV3xowZw8qVK5k4cWLea88//zz9+/fnwQcfTPTpJKnoy8yEF16Ajh3DXWM1a7qgolQEJHwk6Nlnn+Waa64BoE+fPixfvpyuXbvy0UcfJfpUklR8tG8PjzwSHg8bBvfdF2kcSQVQgkqVKpW3cvRxxx1HuXLlAMhx+FdSsuvRA26/PTy+/PKw2aqkyCS8BF144YVcddVVZGdnc8EFF1CxYkXuvfdeOnTokOhTSVLxc/XVYZNVCKXorbeizSMlsYTPCTrjjDPYsGEDJ510Er///e9ZunQpxx57LEOGDEn0qSSp+InFwh5jy5bB889Dp04wbRo0bhx1MinpFOhiiStXrmTfffclFosV1CkKlIslSiowmzaFeUJTp0KNGmExxXr1ok4llQhFYrHEypUrF9sCJEkFqkwZGD8ejjwyjAq1axe+Syo0rhgtSVHZZx94/XWoWxf+/W84+GAYOhR++inqZFJSyHcJGjJkCP/85z9Zvnx5IvNIUnKpVQvefBOaNAnbalx/PRxyCDz2GGRnR51OKtHyPScoJSUl71JXjRo1aNKkCU2bNs37XqtWrYQGjYJzgiQVmpwcePppuO46+Pbb8FqjRmFNoZNPjjabVMwU+C7yzZs3Z+nSpfTu3ZuqVasye/ZsZs2axeeff052djbVqlWjSZMm/POf/8z3HxE1S5CkQrdpE/ztb3DbbbBmTXitffuwwvSRR0abTSomCnxi9IwZM7jlllt4+OGHmTx5MoMGDeKTTz5h3bp1TJ8+nRtvvJHatWvn9/CSlJzKlIErr4Qvv4TLLoNSpcImrI0bwwUXwOLFUSeUSoy4Jkb36tWLf//73xxyyCE0adKEa6+9luzsbJo3b07fvn156KGHEpUzXxYuXEjjxo3zvsqWLcv48eMjzSRJe6RKlbC1xoIF0LVruFw2enTYkPWGG8L8IUlxSdg6QV988QUDBgxg1qxZ3HHHHfTo0SMRh02Y9evXU69ePb755hvKly+/R7/j5TBJRcYHH8DAgTB9enj+u9/BzTfDhRdCWsLXvZWKtUJdJ2jr1q38/PPPnHvuudSuXZvevXuzcuXKRBw6YSZMmMBJJ520xwVIkoqUFi3gvffCbvQHHQQrVkDfvmGe0IQJUHDr3kolVr5L0B133MH5559Pw4YNKVeuHMceeywPPPAAzZo146GHHiIjIyPucNOmTaNjx47UrFmTWCy2y0tZWVlZ1KtXjzJlytC8eXNmzpy5y2M999xznH322XFnkqTIxGLQuTN8+incf3+4ZPb553DGGdC6NfzrX1EnlIqVfJeg6667jg8++IAuXbowf/581q5dy3vvvcfIkSO54IILSE1NjTvchg0baNSoEVlZWbt8/9lnn2XAgAHceOONzJ49m0aNGtGuXTtWrFixw8+tXbuW6dOnc8opp8SdSZIil54eNmH98ku45powmXrqVMjMhPPOg6+/jjqhVCzke07QH/7wB+bOncu6desoX748DRs2pEmTJnlfRxxxREKKUF7QWIyXXnqJTp065b3WvHlzMjMzGTlyJAA5OTnUqVOHfv36cc011+T93OOPP87EiRN54okndnuOn3/+mZ9//jnv+dq1a6lTp45zgiQVbd9+C4MGwRNPhMti6enw17+GNYf23TfqdFKhK/A5QVOnTmXNmjUsXLiQhx9+mOOOO44FCxZwxRVXcNRRR1GhQgWaNWuW38P/ps2bNzNr1izatGmT91pKSgpt2rThgw8+2OFn9/RS2O23305GRkbeV506dRKeW5ISrm5dGDsWZs2Ck06CzZvh7rvhwAPDjvW/+I87SdvFPTG6fv36nHPOOdx1111MnjyZlStX8uWXXzJ27NgdCkqi/fDDD2RnZ1O9evUdXq9evTrLfrEJ4Zo1a5g5cybt2rX7zWNee+21rFmzJu/ru+++S3huSSowRx0FkybBP/8Jhx8Oq1bBgAFw2GHw7LPhNntJefJVgj755BNydvMv0/7770/Xrl0ZOnQoAJ9++ilbt27NX8I4ZWRksHz5ctLT03/zZ0uXLk2lSpV2+JKkYiUWgw4dYO5cePhh2G8/+OorOOecsODi889bhqT/l68SdNRRR/Hjjz/u8c+3aNGCb7fthZMgVatWJTU1dacNXJcvX06NGjUSei5JKnbS0sIaQl98EdYTqlQJ5s2Ds86Chg3DyJAbtCrJ5WuFrdzcXG644QbKlSu3Rz+/efPm/Jxmt9LT02natClTpkzJmyydk5PDlClTuPTSSxN+PkkqlsqXh8GDw91kI0aEVag//TSMDB12WJhQffbZkMAbWaTiIl8lqGXLlixcuHCPf75FixaULVt2r8+zfv16Fi1alPf8q6++Yu7cuVSuXJm6desyYMAAevbsydFHH02zZs2477772LBhA717997rc0lSibbvvnDTTdC/f1hj6N57w5Yc558Pt9wSytA557j6tJJKwrbNKAjvvPMOrVu33un1nj178uijjwIwcuRIhg0bxrJly2jcuDH3338/zZs3T8j53TZDUom1Zk3YrX748DCBGsK+ZIMGhbWGLEMqxvb087tIl6CoWYIklXhr18LIkXDPPbBtu6ODDoLrr4du3SxDKpYKde+wkiYrK4sGDRqQmZkZdRRJKliVKoVFFb/+Gu64A6pWhUWLoHdvOOSQsHP9li1Rp5QKhCNBu+FIkKSks349jBoFw4bB99+H1/bfPxSlHj3CatRSEedIkCRp71WoAFdeGdYWuvtu+N3vwuOLLoKDD4YHHwwrUkslgCVIkrSz8uXhiitCARo+HGrUgG++gT//OcwZGjXK7ThU7MVdgnr27Mm0adMSkUWSVNSUKweXXw7/+U9YZ2i//eC77+AvfwllKCsLNm2KOqWUL3GXoDVr1tCmTRvq16/P0KFDWbx4cSJySZKKkrJlw870//lPuLW+Vi3473/h0kvDRq1jxkSdUNprcZeg8ePHs3jxYvr27cuzzz5LvXr16NChA+PGjWOLdxRIUslSpkwoPosWhVGg2rVhyRLo0wcefzzqdNJeScicoGrVqjFgwAA+/vhjZsyYwUEHHUT37t2pWbMml19+OV988UUiTiNJKirKlAmXxBYtCpfLIEyenjEj2lzSXkjoxOilS5cyadIkJk2aRGpqKqeccgrz5s2jQYMG3HvvvYk8lSSpKChdOtxFdsYZYaL0mWeC0yJUTMRdgrZs2cILL7zAaaedxu9//3uef/55+vfvz5IlS3jssceYPHkyzz33HLfccksi8kqSipqUlHAp7PDDYenSUIQ2bow6lfSb4l4Pfb/99iMnJ4dzzz2XmTNn0rhx451+pnXr1uyzzz7xnqrQZGVlkZWVRXZ2dtRRJKl4qFgRJkyAzEz46KNwaezxxyEWizqZ9KviXjH61ltv5YorrqBcuXI7vJ6bm8t3331H3bp14woYJVeMlqS99PbbcPLJkJ0Nd94JV10VdSIloUJbMfqmm25i/fr1O72+cuVK9t9//3gPL0kqTlq3hvvvD4+vuQZeey3aPNJuxF2Cfm0gaf369ZQpUybew0uSipu+feFPf4LcXDj3XFiwIOpE0i7le07QgAEDAIjFYgwePHiHy2HZ2dnMmDFjl/ODJEklXCwWRoMWLIBp0+D008Ot85UrR51M2kG+S9CcOXOAMBI0b9480n+xs3B6ejqNGjVi4MCB8SeUJBU/6ekwblyYKL1oEZx9Nrz+OqTFfT+OlDBxT4zu3bs3I0aMKJETh50YLUlx+uQTOPZY2LABLrsM7rsv6kRKAoU2MXrMmDEWBEnSrjVsCGPHhscjRsAjj0SbR/qFfI1LDhgwgFtvvZXy5cvnzQ36NcOHD89XMElSCdG5M9x8M9x4Y5g0feihcNxxUaeS8leC5syZk7c56ra5QbsSc5EsSRLAoEEwb16YJ9S5c1hQsRivI6eSIe45QSWZc4IkKYE2bAgjQB9/DI0bw3vvQfnyUadSCVRoc4I2btzITz/9lPf8m2++4b777uPNN9+M99CRycrKokGDBmRmZkYdRZJKjvLl4eWXoVo1mDsXevcOawlJEYl7JKht27Z07tyZP//5z6xevZpDDjmE9PR0fvjhB4YPH07fvn0TlbXQORIkSQXgvffgxBNhyxYYMgSuvz7qRCphCm0kaPbs2ZxwwgkAjBs3jho1avDNN98wduxY7t+2dLokSdscfzw88EB4PGhQGB2SIhB3Cfrpp5+oWLEiAG+++SadO3cmJSWFY445hm+++SbugJKkEujCC+HSS8Pjbt3CpGmpkMVdgg466CDGjx/Pd999x8SJE2nbti0AK1as8BKSJOnXDR8eLoutXw9nnAE//BB1IiWZuEvQ4MGDGThwIPXq1aN58+a0aNECCKNCRx11VNwBJUklVKlS8NxzcMAB8NVX0LVrmCckFZKE3CK/bNkyli5dSqNGjUhJCb1q5syZVKpUiUMPPTTukFFxYrQkFYJPP4VjjgkjQn/5C2RlRZ1Ixdyefn67TtBuWIIkqZC88kq4JJabC6NGwZ//HHUiFWN7+vmdkO18p0yZwpQpU1ixYgU5OTk7vDd69OhEnEKSVJJ17Ai33QbXXQf9+sFhh8Ef/hB1KpVwcc8Juvnmm2nbti1Tpkzhhx9+YNWqVTt8SZK0R665Bs45B7ZuhS5dwjwhqQDFfTlsv/3246677qJ79+6JylRkeDlMkgrZTz9By5YwaxYceSRMnw4VKkSdSsVMoS2WuHnzZo499th4DyNJEpQrB+PHQ/XqYe2gHj3gf6ZZSIkSdwm68MILeeqppxKRpchw7zBJilDt2vDSS5CeHr7ffHPUiVRCxX057LLLLmPs2LE0bNiQhg0bUqpUqR3eHz58eFwBo+TlMEmK0GOPQa9e4fEzz8DZZ0caR8VHod0d9sknn9C4cWMA5s+fv8N7sVgs3sNLkpJVz57wySdhZelu3cLI0JlnRp1KJYjrBO2GI0GSFLHs7FCAnnkG0tLg2Wehc+eoU6mIK7SJ0ZIkFZjUVHj8cTjvvHDr/FlnwbhxUadSCZGQEvTuu+/SrVs3WrRoweLFiwF4/PHHee+99xJxeElSMktLg7Fjw4hQdnZYS+j556NOpRIg7hL0wgsv0K5dO8qWLcucOXP4+eefAVizZg1Dhw6NO6AkSaSmwqOPQvfuoQide264NCbFIe4SNGTIEP7+97/z8MMP73Bn2HHHHcfs2bPjPbwkSUFqKowZE+4Yy84Ol8iefjrqVCrG4i5BCxcupGXLlju9npGRwerVq+M9vCRJ26WmwiOPQJ8+YRHFbt2ghK1Vp8ITdwmqUaMGixYt2un19957jwMOOCDew0uStKOUFHj4YbjwwlCEuneHJ56IOpWKobhL0EUXXcRll13GjBkziMViLFmyhCeffJKBAwfSt2/fRGSUJGlHKSnw4INw0UWhCPXoESZPS3sh7sUSr7nmGnJycjjppJP46aefaNmyJaVLl2bgwIH069cvERklSdpZSgr8/e/bC1GvXpCbGxZZlPZAwhZL3Lx5M4sWLWL9+vU0aNCACiVg118XS5SkYiAnBy69FEaNglgszBnq3TvqVIpQoWybkZOTw6OPPsqLL77I119/TSwWY//99+ePf/wj3bt3d9sMSVLBS0mBrKzt3y+4IBSjCy6IOpmKuHzPCcrNzeX000/nwgsvZPHixRx55JEcfvjhfPPNN/Tq1Ysz3d9FklRYYjH429+gX79wSezCC8PkaWk38j0S9OijjzJt2jSmTJlC69atd3jvrbfeolOnTowdO5YePXrEHbKwZWVlkZWVRXZ2dtRRJEl7KhaDESPCiNCIEXDxxaEQXXxx1MlUROV7TlDbtm058cQTueaaa3b5/tChQ5k6dSoTJ06MK2CUnBMkScVQbi4MGAD33ReejxoFf/5zpJFUuAp8A9VPPvmE9u3b/+r7HTp04OOPP87v4SVJyp9YDIYPD0UIoG9feOCBaDOpSMp3CVq5ciXVq1f/1ferV6/OqlWr8nt4SZLyLxaDu++GgQPD80sugZEjo82kIiffJSg7O5u0tF+fUpSamsrWrVvze3hJkuITi8Fdd8HVV4fn/frB/fdHm0lFSr4nRufm5tKrVy9Kly69y/e37SYvSVJkYjG4/fYwWfr22+Gyy8Lt8/37R51MRUC+S1DPPViRszjeGSZJKmFiMbjttlCEbrsNLr88TJ6+/PKokyli+S5BY8aMSWQOSZIKTiwGt94aitCtt4ZJ0zk5cMUVUSdThOLeQFWSpGIhFoNbboEbbwzPBw6EYcOizaRIWYIkScnlppvCF8BVV8Gdd0aZRhGyBEmSks+NN4ZRIYBrrgmjQt7RnHQsQZKk5HTDDTB0aHh8zz1w8smwYkW0mVSoLEGSpOR17bUwbhxUqADvvANNmsCHH0adSoXEEiRJSm5dusDMmXDoobB4MbRsGfYby9/WmipGLEGSJB12WChCXbrAli3wl79A796wcWPUyVSALEGSJAFUrAjPPx+22khJgcceg+OOg6++ijqZCoglSJKkbWIxuPJKmDQJqlWDOXOgaVN4442ok6kAWIIkSfpfJ54Is2ZBs2awahWcckpYaTonJ+pkSiBL0C5kZWXRoEEDMjMzo44iSYpKnTowbRr86U9hkvTgwXDGGbB6ddTJlCCx3Fynv/+atWvXkpGRwZo1a6hUqVLUcSRJURkzBvr2hZ9/hgMPhBdfhIYNo06lX7Gnn9+OBEmS9Ft694bp0+H3v4cvv4RjjoEnn4w6leJkCZIkaU80aRLmCbVrF26d79YN/vpX2Lw56mTKJ0uQJEl7qkoVeO01GDQoPP/b38Ik6iVLos2lfLEESZK0N1JTw51iL78MlSrB+++H2+jffTfqZNpLliBJkvLj9NPhX/+CI46AZcvCiNCIEW63UYxYgiRJyq/69cOGq+eeC1u3Qv/+cP75sGFD1Mm0ByxBkiTFo3z5cKfYiBGQlgZPPx3uHvvii6iT6TdYgiRJilcsFu4Ue/ttqFED5s+Ho4+GCROiTqbdsARJkpQoxx8Ps2eHjVfXrg0rTN9wg/OEiihLkCRJibTffmFE6K9/Dc+HDIEHHog2k3bJEiRJUqKVKhXmCN11V3h+5ZXOESqCLEGSJBWUK64It85v3Ag9eoQ7yFRkWIIkSSooKSlh89VKlcKt9MOGRZ1Iv2AJkiSpINWtC/ffHx7feCN8/HG0eZTHEiRJUkHr0SPcKbZlC3TvDj//HHUiYQmSJKngxWLw0ENQrRrMmwc33RR1ImEJkiSpcPzud/Dgg+HxXXfB9OnR5pElSJKkQnPmmeHSWE5O+L5+fdSJkpolSJKkwjRiBNSuDV9+CVddFXWapGYJkiSpMO2zT7htHmDUKJg4MdI4ycwSJElSYWvTBi69NDzu0wdWrYo2T5KyBEmSFIU774SDD4YlS6Bfv6jTJCVL0C5kZWXRoEEDMjMzo44iSSqpypWDxx4Lq0o/+SSMGxd1oqQTy83NzY06RFG1du1aMjIyWLNmDZUqVYo6jiSpJBo0CG67DapUgfnzoUaNqBMVe3v6+e1IkCRJURo8GBo3hh9/hIsuAscmCo0lSJKkKKWnw9ix4furr26/c0wFzhIkSVLUjjwSbr01PL7sMvj660jjJAtLkCRJRcEVV8Bxx4VVpHv1CqtKq0BZgiRJKgpSU8PdYuXLw9SpYWVpFShLkCRJRcWBB8I994TH114Ln30WbZ4SzhIkSVJRcvHF0L49/Pxz2GR1y5aoE5VYliBJkoqSWAz+8Q/Yd1+YNQuGDo06UYllCZIkqaipVQuyssLjW2+Ff/0r2jwllCVIkqSi6Jxz4KyzIDs7XBbbuDHqRCWOJUiSpKIoFoMHHgjbaCxYANdfH3WiEscSJElSUVWlSpgfBHDfffDOO1GmKXEsQZIkFWWnngoXXhj2FOvVC9aujTpRiWEJkiSpqBs+HOrVg2++gQEDok5TYliCJEkq6ipWDKtJx2LwyCPwyitRJyoRLEGSJBUHLVtuHwW66CL44Ydo85QAliBJkoqLIUOgQQNYvhz69g3zhJRvliBJkoqLMmVg7FhIS4Nx4+Dpp6NOVKxZgiRJKk6aNoUbbgiPL7kEFi+ONk8xZgmSJKm4ufZayMyE1avDatKbN0edqFiyBEmSVNyUKhUui5UtC2+9BV27WoTywRIkSVJxdOihMH48lC4NEyaEfcYsQnvFEiRJUnHVtm0oQKVLw8svW4T2kiVIkqTirG3bUIC2FaGzz7YI7SFLkCRJxV27dtsvjY0fD+ecA1u2RJ2qyLMESZJUErRvv70IvfRSGBGyCO2WJUiSpJJiWxFKTw9FyBGh3bIESZJUkvyyCL34okVoNyxBkiSVNB06hJGgbUXo3HMtQrtgCZIkqSQ65ZRQgNLT4YUX4LzzLEL/wxIkSVJJdeqp24vQuHFw/vkWoV+wBO1CVlYWDRo0IDMzM+ookiTF59RTw0hQqVLw/POhCG3dGnWqIiGWm5ubG3WIomrt2rVkZGSwZs0aKlWqFHUcSZLy79VXoXPnMBLUtSs89RSkpUWdqkDs6ee3I0GSJCWD005zROh/WIIkSUoWHTtuL0LPPQfduiV1EbIESZKUTDp2DJOkS5WCZ5+F7t2TtghZgiRJSjannx4uiZUqBc88Az16JGURsgRJkpSMzjgjFKG0NHj66aQsQpYgSZKS1RlnhEtj24pQz55JVYQsQZIkJbNfjgg99VQoQtnZUacqFJYgSZKSXadO4W6xJCtCliBJkgRnnhnuFktLgyefhF69SnwRsgRJkqSgc+ftReiJJ0p8EbIESZKk7Tp3DrfNp6aGItSjR4nddNUSJEmSdtSly/YRoaeegj/+ETZtijpVwlmCJEnSzrp0gfHjoUwZmDAh7Ea/fn3UqRLKEiRJknbt1FPh9dehQgV46y1o0wZWrow6VcJYgiRJ0q9r1QqmTIHKlWHGjPB82bKoUyWEJUiSJO1es2YwdSrUqAHz5sEJJ8A330SdKm6WIEmS9NuOOALeew/q1YNFi+D442HhwqhTxcUSJEmS9syBB4YidNhh8N//hhGhuXOjTpVvliBJkrTnatUKl8aaNIHvvw9zhKZPjzpVvliCJEnS3qlWLdwtdvzxsGYNnHwyTJoUdaq9ZgmSJEl7LyMDJk6E9u3hp5/gtNPgpZeiTrVXLEGSJCl/ypWDl1+Grl1h8+awsvRjj0Wdao9ZgiRJUv6lp8PTT0OfPpCTEzZdHTky6lR7xBIkSZLik5oKDz8M/fuH5/36wW23QW5upLF+iyVIkiTFLyUFhg+Hm24KzwcNgquvLtJFyBIkSZISIxaDG2+Ee+8Nz4cNgz//GbKzo831KyxBkiQpsfr3h0ceCaNDDz0E3brBli1Rp9qJJUiSJCVenz7wzDNQqlT4fuaZsHFj1Kl2YAmSJEkFo2vXcAt9mTLw2mtwyimwbl3UqfJYgiRJUsHp0CEsqlixIrzzDpx0Evz4Y9SpAEuQJEkqaC1bwttvQ5Uq8NFHYb+xpUujTmUJkiRJhaBpU5g2DWrWhPnzw75jX30VaSRLkCRJKhwNGsB778EBB8B//gMnnAALFkQWxxIkSZIKz/77w7vvwuGHw+LFcP/9kUWxBEmSpMJVsyZMnQpXXRVpCUqL7MySJCl5VakCd94ZaQRHgiRJUlKyBEmSpKRkCZIkSUnJEiRJkpKSJUiSJCUlS5AkSUpKliBJkpSULEGSJCkpWYIkSVJSsgRJkqSkZAmSJElJyRIkSZKSkiVIkiQlJXeR343c3FwA1q5dG3ESSZK0p7Z9bm/7HP81lqDdWLduHQB16tSJOIkkSdpb69atIyMj41ffj+X+Vk1KYjk5OSxZsoSKFSsSi8V+9efWrl1LnTp1+O6776hUqVIhJlRhy8zM5KOPPoo6RpFUUv7ZFOW/I+pshXn+gj5Xoo+fqOP5eZIYubm5rFu3jpo1a5KS8uszfxwJ2o2UlBRq1669xz9fqVIl/4+2hEtNTfV/419RUv7ZFOW/I+pshXn+gj5Xoo+f6OP5eRK/3Y0AbePEaGkvXHLJJVFHKLJKyj+bovx3RJ2tMM9f0OdK9PGj/t9G+ePlsARYu3YtGRkZrFmzxuYuSco3P08KlyNBCVC6dGluvPFGSpcuHXUUSVIx5udJ4XIkSJIkJSVHgiRJUlKyBEmSpKRkCZIkSUnJEiRJkpKSJUiSJCUlS1AhePXVVznkkEOoX78+//jHP6KOI0kqhs4880z23Xdf/vjHP0YdpcTwFvkCtnXrVho0aMDbb79NRkYGTZs2Zfr06VSpUiXqaJKkYuSdd95h3bp1PPbYY4wbNy7qOCWCI0EFbObMmRx++OHUqlWLChUq0KFDB958882oY0mSiplWrVpRsWLFqGOUKJag3zBt2jQ6duxIzZo1icVijB8/fqefycrKol69epQpU4bmzZszc+bMvPeWLFlCrVq18p7XqlWLxYsXF0Z0SVIREe9niQqGJeg3bNiwgUaNGpGVlbXL95999lkGDBjAjTfeyOzZs2nUqBHt2rVjxYoVhZxUklRU+VlSNFmCfkOHDh0YMmQIZ5555i7fHz58OBdddBG9e/emQYMG/P3vf6dcuXKMHj0agJo1a+4w8rN48WJq1qxZKNklSUVDvJ8lKhiWoDhs3ryZWbNm0aZNm7zXUlJSaNOmDR988AEAzZo1Y/78+SxevJj169fz+uuv065du6giS5KKmD35LFHBSIs6QHH2ww8/kJ2dTfXq1Xd4vXr16nz++ecApKWlcc8999C6dWtycnK46qqrvDNMkpRnTz5LANq0acPHH3/Mhg0bqF27Ns8//zwtWrQo7LgliiWoEJx++umcfvrpUceQJBVjkydPjjpCiePlsDhUrVqV1NRUli9fvsPry5cvp0aNGhGlkiQVJ36WRMcSFIf09HSaNm3KlClT8l7LyclhypQpDlFKkvaInyXR8XLYb1i/fj2LFi3Ke/7VV18xd+5cKleuTN26dRkwYAA9e/bk6KOPplmzZtx3331s2LCB3r17R5haklSU+FlSNLltxm945513aN269U6v9+zZk0cffRSAkSNHMmzYMJYtW0bjxo25//77ad68eSEnlSQVVX6WFE2WIEmSlJScEyRJkpKSJUiSJCUlS5AkSUpKliBJkpSULEGSJCkpWYIkSVJSsgRJkqSkZAmSJElJyRIkqVhp1aoV/fv3z3ter1497rvvvgI9XywWIxaLMXfuXCCs/huLxVi9enWBnfemm27KO29B/n1SMrMESUq4Xr165X2AlypViv3335+rrrqKTZs2JfxcH330ERdffHHCj/tLF110EUuXLuWII44o0PP80sCBA1m6dCm1a9cutHNKycYNVCUViPbt2zNmzBi2bNnCrFmz6NmzJ7FYjDvvvDOh56lWrVpCj7cr5cqVo0aNGgV+nl+qUKECFSpUIDU1tVDPKyUTR4IkFYjSpUtTo0YN6tSpQ6dOnWjTpg2TJk3Ke//HH3/k3HPPpVatWpQrV44jjzySp59+eodjbNiwgR49elChQgX2228/7rnnnp3O88vLYV9//fUOl60AVq9eTSwW45133gFg1apVnH/++VSrVo2yZctSv359xowZE9ff+tNPP9GhQweOO+44Vq9enZfjmWee4dhjj6VMmTIcccQRTJ06dYff+/TTTznttNOoVKkSFStW5IQTTuDLL7+MK4ukPWcJklTg5s+fz/Tp00lPT897bdOmTTRt2pTXXnuN+fPnc/HFF9O9e3dmzpyZ9zNXXnklU6dO5eWXX+bNN9/knXfeYfbs2XFlueGGG/jss894/fXXWbBgAaNGjaJq1ar5Pt7q1as5+eSTycnJYdKkSeyzzz475L/iiiuYM2cOLVq0oGPHjvz4448ALF68mJYtW1K6dGneeustZs2aRZ8+fdi6dWtcf5+kPeflMEkF4tVXX6VChQps3bqVn3/+mZSUFEaOHJn3fq1atRg4cGDe8379+jFx4kSee+45mjVrxvr163nkkUd44oknOOmkkwB47LHH4p4j8+2333LUUUdx9NFHA2EkKb+WLVvG2WefTf369Xnqqad2KHkAl156KV26dAFg1KhRvPHGGzzyyCNcddVVZGVlkZGRwTPPPEOpUqUAOPjgg/OdRdLeswRJKhCtW7dm1KhRbNiwgXvvvZe0tLS8QgCQnZ3N0KFDee6551i8eDGbN2/m559/ply5cgB8+eWXbN68mebNm+f9TuXKlTnkkEPiytW3b1+6dOnC7Nmzadu2LZ06deLYY4/N17FOPvlkmjVrxrPPPrvLuTstWrTIe5yWlsbRRx/NggULAJg7dy4nnHBCXgGSVPi8HCapQJQvX56DDjqIRo0aMXr0aGbMmMEjjzyS9/6wYcMYMWIEV199NW+//TZz586lXbt2bN68Od/nTEkJ/y8tNzc377UtW7bs8DMdOnTgm2++4fLLL2fJkiWcdNJJO4xI7Y1TTz2VadOm8dlnn+3175YtWzZf55SUOJYgSQUuJSWF6667jkGDBrFx40YA3n//fc444wy6detGo0aNOOCAA/j3v/+d9zsHHnggpUqVYsaMGXmvrVq1aoef+V/b7hRbunRp3mu/nCT9y5/r2bMnTzzxBPfddx8PPfRQvv6uO+64g549e3LSSSftsgh9+OGHeY+3bt3KrFmzOOywwwBo2LAh77777k4lTVLhsQRJKhRdu3YlNTWVrKwsAOrXr8+kSZOYPn06CxYs4E9/+hPLly/P+/kKFSpwwQUXcOWVV/LWW28xf/58evXqlTfasytly5blmGOO4Y477mDBggVMnTqVQYMG7fAzgwcP5uWXX2bRokV8+umnvPrqq3nFJD/uvvtuzj//fE488UQ+//zzHd7LysripZde4vPPP+eSSy5h1apV9OnTBwjzhdauXcs555zDv/71L7744gsef/xxFi5cmO8skvaOJUhSoUhLS+PSSy/lrrvuYsOGDQwaNIgmTZrQrl07WrVqRY0aNejUqdMOvzNs2DBOOOEEOnbsSJs2bTj++ONp2rTpbs8zevRotm7dStOmTenfvz9DhgzZ4f309HSuvfZaGjZsSMuWLUlNTeWZZ56J62+79957OeusszjxxBN3GKm64447uOOOO2jUqBHvvfceEyZMyLsTrUqVKrz11lusX7+eP/zhDzRt2pSHH37YOUJSIYrl/vLiuSRpB61ataJx48Z7tXXF119/zf7778+cOXNo3LhxXOevV68e/fv332GrEEmJ4UiQJP2GBx54gAoVKjBv3rxCO+fQoUOpUKEC3377baGdU0o2jgRJ0m4sXrw4bzJ33bp1d1oLaFcSMRK0cuVKVq5cCYSJ3BkZGfk6jqRfZwmSJElJycthkiQpKVmCJElSUrIESZKkpGQJkiRJSckSJEmSkpIlSJIkJSVLkCRJSkqWIEmSlJQsQZIkKSn9HyHgjM5HOcMxAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# now, we can make this prettier according to our use-case:\n", "import matplotlib.pyplot as plt\n", "fig, ax = plt.subplots()\n", "ax.plot(bins, density, c='r')\n", "ax.set(xscale=\"log\", yscale=\"log\", xlabel=\"Radius [kpc]\", ylabel=r\"Density [$M_\\odot / kpc^3$]\")" ] }, { "cell_type": "markdown", "id": "e51ae5b7", "metadata": {}, "source": [ "As another example, here's a surface density plot of the particles we selected:" ] }, { "cell_type": "code", "execution_count": 15, "id": "06b3bcdd", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbUAAAGiCAYAAABgeVj+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0XklEQVR4nO3dfXBV1b3/8U8i4ZAISYgkBFoepWKtikhrGh86UBkeruOFynXsw7Vg/eGIQUeCVdJWFLwalFa9tQi991eBTkexTke9Ol6vSgHHS9DKNdenQgXBgJKgIImgngTYvz/89bj3Osne2TkP+2Txfs2cmezsfdb+nn2CX/f67rVWnuM4jgAAsEB+1AEAAJAuJDUAgDVIagAAa5DUAADWIKkBAKxBUgMAWIOkBgCwBkkNAGANkhoAwBokNQCANTKa1Orr6/Wtb31LAwYMUEVFhWbOnKnt27d7jvn8889VU1OjU045Rf3799esWbPU0tKSybAAAJbKaFLbtGmTampqtGXLFj3//PPq6OjQlClTdOTIkcQxCxYs0FNPPaXHHntMmzZt0gcffKDLLrssk2EBACyVl80JjT/88ENVVFRo06ZN+s53vqPW1laVl5fr4Ycf1j/90z9JkrZt26avf/3ramho0Le//e1shQYAsECfbJ6stbVVklRWViZJ2rp1qzo6OjR58uTEMaeffrqGDx/eZVKLx+OKx+OJ7ePHj+vgwYM65ZRTlJeXl+FPAABIN8dx9Mknn2jo0KHKz0+tAzFrSe348eO68cYbdcEFF+jMM8+UJDU3N6tv374qLS31HDt48GA1Nzd32k59fb2WLFmS6XABAFm2Z88effWrX02pjawltZqaGr355pt66aWXUmqnrq5OtbW1ie3W1lYNHz5cF+of1EcFqYYJAMiyo+rQS3pGAwYMSLmtrCS1+fPn6+mnn9aLL77oycKVlZVqb2/XoUOHPHdrLS0tqqys7LStWCymWCyW9Ps+KlCfPJIaAPQ6///JjnSUkDL69KPjOJo/f74ef/xx/fnPf9aoUaM8+ydMmKCCggKtX78+8bvt27erqalJ1dXVmQwNAGChjN6p1dTU6OGHH9aTTz6pAQMGJOpkJSUlKiwsVElJia6++mrV1taqrKxMxcXFuv7661VdXc2TjwCA0DKa1FauXClJmjhxouf3q1ev1pw5cyRJ9913n/Lz8zVr1izF43FNnTpVDz74YCbDAgBYKqvj1DKhra1NJSUlmqgZ1NQAoBc66nRoo55Ua2uriouLU2qLuR8BANYgqQEArEFSAwBYg6QGALAGSQ0AYA2SGgDAGiQ1AIA1SGoAAGuQ1AAA1iCpAQCsQVIDAFiDpAYAsAZJDQBgDZIaAMAaJDUAgDVIagAAa5DUAADWIKkBAKxBUgMAWIOkBgCwBkkNAGANkhoAwBokNQCANUhqAABrkNQAANYgqQEArEFSAwBYg6QGALAGSQ0AYA2SGgDAGiQ1AIA1SGoAAGuQ1AAA1iCpAQCsQVIDAFgjo0ntxRdf1KWXXqqhQ4cqLy9PTzzxhGf/nDlzlJeX53lNmzYtkyEBACyW0aR25MgRjRs3TitWrOjymGnTpmnfvn2J1yOPPJLJkAAAFuuTycanT5+u6dOn+x4Ti8VUWVmZyTAASOozxPvv7Oi+5ogiATIn8praxo0bVVFRobFjx2revHk6cOCA7/HxeFxtbW2eFwAAUsRJbdq0afr973+v9evX6+6779amTZs0ffp0HTt2rMv31NfXq6SkJPEaNmxYFiMGAOSyPMdxnKycKC9Pjz/+uGbOnNnlMe+++65OPfVUvfDCC7r44os7PSYejysejye229raNGzYME3UDPXJK0h32IA16H5ErjrqdGijnlRra6uKi4tTaiujNbWwRo8erUGDBmnHjh1dJrVYLKZYLJblyIDcYCYmt6AkRRLDiSDymprb3r17deDAAQ0ZMiTqUAAAvVBG79QOHz6sHTt2JLZ37dqlxsZGlZWVqaysTEuWLNGsWbNUWVmpnTt36uabb9aYMWM0derUTIYFALBURpPaq6++qkmTJiW2a2trJUmzZ8/WypUr9frrr2vt2rU6dOiQhg4dqilTpuiOO+6gexEnrKC6F12IgL+MJrWJEyfK7zmU//qv/8rk6QEAJ5icqqkBAJAKkhoAwBo59Ug/YKMw48OCamZ+bfk97t9Z2+7jU6nVMf4NuYQ7NQCANUhqAABrkNQAANagpgZkWDprTKnU41I9vrvthKn7UX9DunGnBgCwBkkNAGANkhoAwBrU1IAMy1YdKVfqVdTJECXu1AAA1iCpAQCsQfcjThhRdc+FeeQ9SP45Z3i2jze+3eO43G2l0k4YdE0i07hTAwBYg6QGALAGSQ0AYA1qajhhZKueE1QzCxNHUls+ta/Qy9ZkqY7mGwM1NqQZd2oAAGuQ1AAA1iCpAQCsQU0NGRGmdpLJOou77aB2/cZtmePD8lsOerYzVRvK1rI1plS+E7/3pvJ5qMehO7hTAwBYg6QGALAG3Y/IiDBdQ5nqbgzbtt9UUUHdjX7dnOmc2ipbcrFrLxdjQu7hTg0AYA2SGgDAGiQ1AIA1qKkh54V5LD9dj4wfH1zmbTegDuZ33kzW0MIsH2PWG92fMWyMft+J37XomPJNz3bhG3u7/V6/GDp7b5i/G9iDOzUAgDVIagAAa5DUAADWoKaGHsnUNEqdce8Pem8qY8I85zH2He92K9nl/nxhr6vfe02pLJ/jVvDcq952etRK92KgjnZi4k4NAGANkhoAwBoZTWovvviiLr30Ug0dOlR5eXl64oknPPsdx9HixYs1ZMgQFRYWavLkyXrnnXcyGRIAwGIZrakdOXJE48aN009+8hNddtllSfvvuece/frXv9batWs1atQo3XrrrZo6darefvtt9evXL5OhIUWp1CvCvjfU8jHGHI2ZqoXlytI6Ydq1fRwXS9NAynBSmz59uqZPn97pPsdxdP/99+sXv/iFZsyYIUn6/e9/r8GDB+uJJ57Q97///UyGBgCwUGQ1tV27dqm5uVmTJ09O/K6kpERVVVVqaGjo8n3xeFxtbW2eFwAAUoSP9Dc3f9E1MHjwYM/vBw8enNjXmfr6ei1ZsiSjsSH9UunqSqUbqafTSPWWx8WzdV3TOYSjp+0ECdMFnEpbyG297unHuro6tba2Jl579uyJOiQAQI6ILKlVVn7xf04tLS2e37e0tCT2dSYWi6m4uNjzAgBAijCpjRo1SpWVlVq/fn3id21tbXr55ZdVXV0dVVgAgF4sozW1w4cPa8eOHYntXbt2qbGxUWVlZRo+fLhuvPFG/cu//Iu+9rWvJR7pHzp0qGbOnJnJsBABv6muwrzXFFjrCVEbSdeyNalMKRVmOZUg5jIv5hRVYYT5DvyW7QmKP52P5WerlofcktGk9uqrr2rSpEmJ7draWknS7NmztWbNGt188806cuSIrrnmGh06dEgXXnihnn32WcaoAQB6JKNJbeLEiXIcp8v9eXl5Wrp0qZYuXZrJMAAAJ4he9/QjAABdYekZJORKPSNM3ShT9apUYkzlOgbV0HyX2jH2Fb6x19t4Gqf2cp/3qDn+z6fdTE7l5VfbCzNGsTvnQu7iTg0AYA2SGgDAGiQ1AIA1qKkhIZ01tEwtTZPOelUYqSw146l7KXl5HL+2AsepudrKDxgvFlRX6ul5zeV9zM8b5rxBcYQ51v0fN/PapDKWELmNOzUAgDVIagAAa9D9iB4JO9VVKm27u4bCPtadru5J30fp5e1S9OsGk5K7BcNM5eXXpWjGFEbQdQtzHeMVRZ7tghTi8BP0nSQNNcAJgTs1AIA1SGoAAGuQ1AAA1qCmhm4LM2VRKsJMlZTOoQT755+f+Hng39o9+8xlW8zH1M3H2v1iSqqxdTvC5OEAnkfVjX1J781SrdJviZuw31eY79Nv6EAmp8Fiiq3cwp0aAMAaJDUAgDVIagAAa1BTO8GkUhdzH5/NOkKY8WOmMJ936J/e7fo8Pu2abfstgSIpsPYVhrtts972wazRnu2K32z2bLtriOZnz9Q0UqnURIPGsKUy3VimpnVD9nGnBgCwBkkNAGANkhoAwBrU1CwTZg6/bC7jkkotzz2nX9CyLal8ps/O+mriZ3OsVdK8ij51saBxaeb+jinfTPxc+MZe32OT5jt01dTMfUP/ZJzYuDaVLx1Sd/nVs8KOcQvDfd3TOZdjOmuxyC3cqQEArEFSAwBYg+5Hy2Rr+p9MPtLvt3JymCmlJG+c5qP15hIpnq4/n5WepW6sDO3D/Hyx/Z922a67a9I8VgpYesaI+eCkkZ7t4oe3fNmOEWM6v1+/YQdh2g3bjZnK32umVl9H5nGnBgCwBkkNAGANkhoAwBp5juM4UQeRira2NpWUlGiiZqhPXpiF42Hyq2WlUzprEplqK6hm5n78X/LW48LW21KpOblrbma9zWwrabouF7O+GDSkwf23kcr0Velkxmh+fj+ZijHo75OhA1846nRoo55Ua2uriouLU2qLOzUAgDVIagAAa5DUAADWoKZmgUz1y0c1HieVZV2aLyz1bPsuJxNiLJJ5HrPe6D7er94WFEcq9amw47jcnyls/dTvXH7jAc1aXZjzhJ2Oi7FmvQc1NQAAOkFSAwBYg6QGALAGcz/2AmGWk0mndNZvcrG+Eea8Zs3JHBPlXhYl5jMeLOg8gXM/+rUdMC7Nb9xW2O/Ab8ybn7A1w3TNN5rKe83vwKwLUp/LLZHfqd1+++3Ky8vzvE4//fSowwIA9EI5caf2jW98Qy+88EJiu0+fnAgLANDL5ET26NOnjyoru9edFY/HFY/HE9ttbW2ZCitn5Ep3RlRT+vh1WZldahW/MR5NT1PMSV1MPo/AB03PZHZn+T3y3ydgJWzfR+uNGJOWl/GJMYj7MwZd18IQj/8HrRLux7f72OwuDbFietAwhFz594kvRN79KEnvvPOOhg4dqtGjR+tHP/qRmpqaujy2vr5eJSUlidewYcOyGCkAIJdFntSqqqq0Zs0aPfvss1q5cqV27dqliy66SJ988kmnx9fV1am1tTXx2rNnT5YjBgDkqsi7H6dPn574+eyzz1ZVVZVGjBihP/7xj7r66quTjo/FYorFYtkMEQDQS0Se1EylpaU67bTTtGPHjqhD6TXcj5f7Td8kpdb/n63Hq8MMFzCnpCow2vLUUtJUn5HCff6Dk0Z6tosf3uI9IETNyfwMZttuZRu63NV52z7nSYX72phDIfyODdoXNJWZp+4XMA1YpoYH8Lh/9kXe/Wg6fPiwdu7cqSFDhkQdCgCgl4k8qd10003atGmTdu/erc2bN+t73/ueTjrpJP3gBz+IOjQAQC8Teffj3r179YMf/EAHDhxQeXm5LrzwQm3ZskXl5eVRhwYA6GUiT2rr1q2LOoRez2/ZkDB9+Ga9w68+l67aXHf2u+Myx4CZY4j8ajZhpzvy41frMuMv/Mhb9ws7xVh3mbW6o8b+MOcNOtZdyzSn7vKbUsz8/sLUU836qczxY2Y9tcuWwtV1g6al8/t3k62pvPClyLsfAQBIF5IaAMAaJDUAgDXyHMdxog4iFW1tbSopKdFEzVCfvIKow0mLsOO4erpETNg+e7/xcEH1uHQJu3SJX40tqL7j1475XneNLWjuR7MeF68oSvxszgNp1pE+G+Qtg5dt2J34OZXvs/nCUs++it9s9mz7zVdpxhjb/6ln2/23EDROLWicpZt5Xr85NINQv4rWUadDG/WkWltbVVxcnFJb3KkBAKxBUgMAWCPyR/qRLJ1dIalMQWVydw0FLcXiuyRKwCPgft2rYT/PoTO+7MowH3nPN94bZshCUhyun4O6DM04Cn3Oaz4uH9vvjcPT7Rmw0rXftFgD/9be5T5J+vi0vkYcXbdlfp/urkuza9KMMW50c+a7jjfbTWWF8TDCdnkjWtypAQCsQVIDAFiDpAYAsAY1tV4ok1NU9bQd81FtvxpbUq0noG33e826kHkeMw7PI+/y5267zVjSxd2OlFyPS5rCyaXwI++Z/ZZMMf9BfuZ63F9Krs/5ncdcddCvxpZU6zJirHzpUJfnDZpuzF0XCxpGUdjlWZJjMr97v2mxUpHOujQyjzs1AIA1SGoAAGuQ1AAA1qCm1guka1os873pXOoibJ3MzYzDrE/1cU1/FDiNkt/YLHNpEp/xcmUtXde9OovDXQsKrDca+/fPGp34eeDf/GtopW+3edt2fV4zxjDj1ILEjdqeuwYX9P25FYacQsxvurVsLdXCkjC9C3dqAABrkNQAANag+7EXCNPdkaljpXAz/vt12QR15xSY3YTu7ruAVYhN7j9wcwom81F0zzkDusHMWes73Bs+3W+SFDO6I91TVJndjeZj+n7cU4JJUqmx3+yOPOgatmCex+wm9Jt5X8bnCTNbfmD3sevnoGEj6RRmyrQwK1SkUkqg27N7uFMDAFiDpAYAsAZJDQBgDVa+tkC6VrdOZfqfoJpaKitDu49PpQ5ovrfth9/2bJuPy/vx+wxm/O9dWurZjpd5BzyctvbL85p1sc/LvP/f6TddlfnYvblcjPle9/FmfTFohep0CbOCut/fVGdtuVGfym2sfA0AQCdIagAAa5DUAADWYJxaL+RXW/CrK5jvNesK6awzJLXl2g5aIiTfXF7GZ4ybyZyiyb3sSYexL5Ua2kGfpWnMY0c85d/WB65psvod9F6d+EDve826mR+/Gprkrbn1GxSuvuhuK6h25xZUQzNjdFfJk+q2ZtvGtmcsoc+YRCl9YzDN/anUtKn79Qx3agAAa5DUAADWIKkBAKxBTa0X8qtXhVl6JkiY+RuDhFoux6g55Yc4l9+8gwVGPcecC9I9v6E5XqzYeG/p2966n7vGZtajzHFqo/6v9/O562YD/+adg/HzMm+9yuSeK9KcvzFo3Jp7zklzbsegZWvc17ngOeO7NepknuVx5M/8/ro/82Uyv7k9TemcM7WntTBqaOnBnRoAwBokNQCANeh+tIDfI8TpnPqqu/tS5de2+Qi42U3m9/mDpn5yd9e5H9GXjOVvOjlvmWvbHFZQ+XK7Z9scDhD7+MufzS5C8xH/5irj8XlX22Z3m/ndf3zaaM+2u6svacqpEKtkJ3VTm4/tu/YnDUcx2srUUixBq3P7dVWGPW9PH+lHenCnBgCwBkkNAGCNnEhqK1as0MiRI9WvXz9VVVXplVdeiTokAEAvFHlN7dFHH1Vtba1WrVqlqqoq3X///Zo6daq2b9+uioqKqMPrFcLUvsL096drSZuw5zW5a2FBy9b4xWy+13zE36/mFvSIe/OFpYmfh/7pXc8+s35jPvK//cbCxM8DSr2P1l84fJtn+4WHqj3bE5Zt/XLfaed79rWd6q1YxYxL51c3M4cDuIcOSFKpux3jOoYZCpLKtG6p1NjMoQNmzdQdVyr/TpB9kd+p3XvvvZo7d66uuuoqnXHGGVq1apWKior00EMPdXp8PB5XW1ub5wUAgBRxUmtvb9fWrVs1efLkxO/y8/M1efJkNTQ0dPqe+vp6lZSUJF7Dhg3LVrgAgBwXaVL76KOPdOzYMQ0ePNjz+8GDB6u5ufNb/rq6OrW2tiZee/bsyUaoAIBeIPKaWlixWEyxWCzqMLIqneNzQk1X5SPseDG/ukrQ5/HUXUKet+2HXy6pUvzwFs8+v+m5zHqTOW3W52Xe/x9019H8lqWRpG03e/cPfs7drnetmecmn+7Z/t+fPejZHn/XdYmfK36z2bOv0rhWf5vt/Qxu5uczp9zyW4qmw5hurI9Rr3L/rZhj2HJFUG3PT7rGorH0THpEeqc2aNAgnXTSSWppafH8vqWlRZWVFF8BAOFEmtT69u2rCRMmaP369YnfHT9+XOvXr1d1dbXPOwEASBZ592Ntba1mz56tb37zmzrvvPN0//3368iRI7rqqquiDg0A0MtEntSuuOIKffjhh1q8eLGam5t1zjnn6Nlnn016eORElsr4nDDj1MKMxzFrEGGWhzHrVUF/hO64wo5Tc9fRgmqGfuO2zHpcqVGv+mDWl/MqmuPUzBpbxV+8bbvrVeb4sPdOLfFsV//xWs92P9fsifvne8epmXHEDpZ6tt11NLPuZzLH2rmXqjHHfJnHuhUG/O2mUlcy67x+dbJ0ntcvjjC1urDj4ai5dS7ypCZJ8+fP1/z586MOAwDQy0U++BoAgHTJiTs1hBNVt0O6lqIxu2SCVkP2HBuwJIrZPenpkgrouvRrJyjGypcOJX42u9/Mrr1d/8dYAuajL7sczRWoi3eWGmfyRuIeWmAuU2NeK3MJHLeg7jgzLnc3qTnAJmlKLZ/hAEFdap4ub3N5HKOtpO8+hSWZ/Lrpg7o5w6z0jfTjTg0AYA2SGgDAGiQ1AIA1qKnloKD+/jD1gTDTVSXVLAKWFPHjV/sy2zWnWSp47tUuzxtUR0mq7bm2D7mmzJI6mTbLp52gacHczBqT+Uh/yc6uKy3mdFUD/9Z1HUySSt/+tMt97uVwOmvL8yh+wHdrft5C17bfI/yS9/sOqjGZ19k9rZb5HyvzO0qlfhVq6EDIeqtbKo/l8wh/93CnBgCwBkkNAGANkhoAwBrU1HJQ2DE1fu8NI2hKn1DLx4SIw6yhmXUVz5IwAWPczPqcZzonYzmVVJg1Q3N6KzfzvOaYL3eNJhbQrjkllXt6LvdYuc62zRpbbP+X5wpbP3V/34XGvs+MGmKYaaPCTG0VJOnvyN2WObYshSm2wu7vLqbF6hnu1AAA1iCpAQCsQVIDAFiDmlovkKllMlJZ6iKV/v2gcWlBy8u4mTHm+9SrCs3xRT61O7Pdz4zallkXc2+btSuztmWORStzndf87B8bbX02aKRn21xexs0cH+dXY6to9L43aMyiH3P8n7vuGXYMpvtvxawnBr3XL+akv+0wx6ZxKSg/1NB6hjs1AIA1SGoAAGvQ/dgLhemWyNS0PGGmjZK83VlBj/CbXUHuzxD0efKNbjPfz2u+19W2OfWTGbN8Pv/QPxndnAHL5bhjNLtmK36z2bNt7vc7j7nki9nt6e66NLtigx7x95wrRHdx0PRjHT4rbAfJ1CPw6ewG5DH9zONODQBgDZIaAMAaJDUAgDWoqVnAb+n5TPXZB01XlcRnii2zJmO25VdzKjBrNCksC+Jp16ih+U2/FcSc6sqsdblreX3MabDmn+/ZNh/Ld9f+zEfeze++uNEb13Gf6avC1i49+3zqc0HffYE5fVWIOme2ppNLqmv6xJHKeai/9Qx3agAAa5DUAADWIKkBAKxBTa0XCOpbD7UkTBrP68evrpJKTGHrKH7LjyTVBd11ooAYQ9UUjVqQ35g+s27kNw2WlDz1l+c8AWPPPHH4TBkWJKje5o7DbwxiUNvmsjxmzH5TXaWTWU/1qwGnc5woNbbu4U4NAGANkhoAwBokNQCANaip9QK5uFxF4JggoyaT73N8mCVwQo/Dc+03a0x+y48ELmnjMxdiUs3MqAOa8yya49jcPpg12rNtzgWZ7zfWzGjLt/bTZQSdc39Gs92ksYKua5PKMi5hlr8x2wpaSsfv7znMkjZBbaWCGlr3cKcGALAGSQ0AYA26H9Ft6Zr+xxT06HmmltpJmnYpxHnM7ixPd5yxz/xHdtBYAsacNsutMqDb093Ne+iH3/bsKtuw2xtHGqeR8kzPZezzuzamVLrq0vmIu997w57H75F+ZB53agAAa5DUAADWIKkBAKwRaU1t5MiReu+99zy/q6+v16JFiyKKqHfKxaUu/Nrye5RektqM2lDxw1vSEqNZqys02+p2S/51oqTpqIzzmLUu3+EAZts+j6aXbfA/1u9Re3NfYC3INUzBHKIQZgqxoHqqn6DhAH6fz2+ISdB5kNsif1Bk6dKlmjt3bmJ7wIABEUYDAOjNIk9qAwYMUGUlTwgBAFIXeU1t2bJlOuWUUzR+/HgtX75cR4/6dwLF43G1tbV5XgAASBHfqd1www0699xzVVZWps2bN6uurk779u3Tvffe2+V76uvrtWTJkixGaa9MjqHxq0MkLVVi1FXcNbQw7XZnv9+xftMqmfzqVebnCbPcijl2rvCNvd7zmvU697nMupGxnfR53e0GfV6f8ybV6nzaSfrsId4bJOm7d0+xFbIt9+cL+x9JdxypTAGHnkn7ndqiRYuUl5fn+9q2bZskqba2VhMnTtTZZ5+ta6+9Vr/61a/0wAMPKB6Pd9l+XV2dWltbE689e/ak+yMAAHqptN+pLVy4UHPmzPE9ZvTo0Z3+vqqqSkePHtXu3bs1duzYTo+JxWKKxWKphgkAsFDak1p5ebnKy8t79N7Gxkbl5+eroqIizVHZradT/AQ9Eu3XVphuFcnbvZX0yHcKKxan8ki4Xzdh0urNAY/ae66HcW06pnzTs93Hp0uxwGj3M+O95lRe7rjCfp9dxdAZv+tsdu2Zn9e9WnlQjJkaVhK2nTDXNZUptuhyTL/IamoNDQ16+eWXNWnSJA0YMEANDQ1asGCB/vmf/1kDBw6MKiwAQC8WWVKLxWJat26dbr/9dsXjcY0aNUoLFixQbW1tVCEBAHq5yJLaueeeqy1bun7KDQCAsPIcx3GiDiIVbW1tKikp0UTNUJ88syIB25m1nqT9Paw5BdU6zPP61dyCHtMPw6+toKESJr/lcoKGA7ilUgdLZw3NNifStTnqdGijnlRra6uKi4uD3+Aj8sHXAACkC0kNAGANkhoAwBqRT2iM9ErnGJpsSSnmgPFVftMjpfJ5/ZZXMWM0a2hhPp/Jrx5n1tDiFUWe7dj+T7t+b4gamuQ/BizM+Lh01jlTEWqppCz9O7G5hpZJ3KkBAKxBUgMAWIOkBgCwBjU1y6Rr6ZXu7O/pecPGka7zhPl86fw8YePwEzSezM0956Ik5Yc4b5g5NcNcVzN+8z9AqVz3VOYBDbNUEnIbd2oAAGuQ1AAA1qD7ET2SyW7ATJ3LbMdvSZigx9IPThrp2S59uy3xc6hlakLsk/y768x95nCGUEMJzFW0fb4zs9vP7z8q5jCDwoBrFUaY7kam9rIXd2oAAGuQ1AAA1iCpAQCsQU0NCZmsbaXrPEG1LvORcb9ptMxpo8Kct/Cjo97zus7jNzVXWGbdT+Zj+q56VtCUYWGmrzKvo1mv8rQVYoqtwm5HkHyeTA6VCIMaWm7jTg0AYA2SGgDAGiQ1AIA1qKlZJp1TQYWpUZhjlcIuZdJdge2E+Lxmncj9GYLiL0jjuCa/GJOWi/GZCiroPKHGYhnbYabY8osj7N9BKvXVnraL3o07NQCANUhqAABrkNQAANagpmaZdI41C7UkijlGKsR8gOmMI5U6i/szmDUkk1lD9OwLqMeZ7z3qVxfzGx8mea5z0Di1MHJxLsRUlkbCiYM7NQCANUhqAABr0P14Astkt5HZHZnK1FF+Kyeb/LrgwnzeoGvh18Vodl0GdSl6psIypsEKisv9DzioyzRTK4yn8+8mU6tXmzL5t5+uFdTRM9ypAQCsQVIDAFiDpAYAsAY1tRNYVP39YesmnjhDxtzTacGytQyPJBW46mhhaz3u/WEfec/Feo/f1GXmEjepyORnz8XreiLhTg0AYA2SGgDAGiQ1AIA1qKkh64JqaKnUfrI1VVI6l1dJ5bxdxRC2rbBTTvnV8rr7vm6dJ8TUZYCUwTu1O++8U+eff76KiopUWlra6TFNTU265JJLVFRUpIqKCv30pz/V0aNHMxUSAMByGbtTa29v1+WXX67q6mr97ne/S9p/7NgxXXLJJaqsrNTmzZu1b98+/fjHP1ZBQYHuuuuuTIUFALBYxu7UlixZogULFuiss87qdP9zzz2nt99+W3/4wx90zjnnaPr06brjjju0YsUKtbe3ZyosAIDFIqupNTQ06KyzztLgwYMTv5s6darmzZunt956S+PHj+/0ffF4XPF4PLHd1taW8ViRXkE1M7/akPlecy5Iv7kfg8bHhRnzFaZ+lcp8hqZMzW8YpoaWSkypvDeb479SqRMiWpE9/djc3OxJaJIS283NXf+R1NfXq6SkJPEaNmxYRuMEAPQeoZLaokWLlJeX5/vatm1bpmKVJNXV1am1tTXx2rNnT0bPBwDoPUJ1Py5cuFBz5szxPWb06NHdaquyslKvvPKK53ctLS2JfV2JxWKKxWLdOgeiE+aR91Qe4T90RrFnu7ix62PDDCUIXHrGp0vRb1XsoPMG7fO7dmG7SP2moErnlGJdvS+T5+nsXH7n7Q1TiKF7QiW18vJylZeXp+XE1dXVuvPOO7V//35VVFRIkp5//nkVFxfrjDPC/UcBAAApgw+KNDU16eDBg2pqatKxY8fU2NgoSRozZoz69++vKVOm6IwzztCVV16pe+65R83NzfrFL36hmpoa7sQAAD2SsaS2ePFirV27NrH996cZN2zYoIkTJ+qkk07S008/rXnz5qm6ulonn3yyZs+eraVLl2YqJACA5fIcx3GiDiIVbW1tKikp0UTNUJ+8gqjDQRfSNb1T2OVV0nWeVNoKajvMo/apnDfM0IJsDQcAJOmo06GNelKtra0qLi4OfoMPJjQGAFiDpAYAsAZJDQBgDZaeQVakq64S1E62zuMn7PI3fmOxzGnAZOx318mCxuGFmZ4rlfFj2Vr+B+gMd2oAAGuQ1AAA1iCpAQCsQU0NvUo6l3HJFLPGlErM6ayT+Y09CzserqdL0TDHIjKNOzUAgDVIagAAa9D9iF4lzPIxUm50b4WN2S1b8WdrCEMufB+wG3dqAABrkNQAANYgqQEArEFNDZFLZx2sN9Rsoqr7Zeo8veGa48TBnRoAwBokNQCANUhqAABrUFND5HKlJhNmOZV01v2ydd5cHMMHpBt3agAAa5DUAADWoPsRVumNXWzpijHos/eGawGkijs1AIA1SGoAAGuQ1AAA1qCmBqtkaoqtoHb9HstPRZhH/KmZAdypAQAsQlIDAFiDpAYAsAY1NaATQWO+bFs+xtQbx/sBEndqAACLkNQAANYgqQEArEFNDehEUA0pnTWmTI1xS0WuxAGExZ0aAMAaGUtqd955p84//3wVFRWptLS002Py8vKSXuvWrctUSAAAy2Ws+7G9vV2XX365qqur9bvf/a7L41avXq1p06YltrtKgICt6OoD0idjSW3JkiWSpDVr1vgeV1paqsrKSt9jAADojshrajU1NRo0aJDOO+88PfTQQ3Icx/f4eDyutrY2zwsAACnipx+XLl2q7373uyoqKtJzzz2n6667TocPH9YNN9zQ5Xvq6+sTd4EAALiFulNbtGhRpw93uF/btm3rdnu33nqrLrjgAo0fP1633HKLbr75Zi1fvtz3PXV1dWptbU289uzZE+YjAAAsFupObeHChZozZ47vMaNHj+5xMFVVVbrjjjsUj8cVi8U6PSYWi3W5DwBwYguV1MrLy1VeXp6pWNTY2KiBAweStAAAPZKxmlpTU5MOHjyopqYmHTt2TI2NjZKkMWPGqH///nrqqafU0tKib3/72+rXr5+ef/553XXXXbrpppsyFRIAwHIZS2qLFy/W2rVrE9vjx4+XJG3YsEETJ05UQUGBVqxYoQULFshxHI0ZM0b33nuv5s6dm6mQAACWy3OCnqHPcW1tbSopKdFEzVCfvIKowwEAhHTU6dBGPanW1lYVFxen1Fbk49QAAEgXkhoAwBokNQCANUhqAABrkNQAANYgqQEArEFSAwBYg6QGALAGSQ0AYA2SGgDAGiQ1AIA1SGoAAGuQ1AAA1iCpAQCsQVIDAFiDpAYAsAZJDQBgDZIaAMAaJDUAgDVIagAAa5DUAADWIKkBAKxBUgMAWIOkBgCwBkkNAGANkhoAwBokNQCANUhqAABrkNQAANYgqQEArEFSAwBYg6QGALAGSQ0AYA2SGgDAGiQ1AIA1SGoAAGtkLKnt3r1bV199tUaNGqXCwkKdeuqpuu2229Te3u457vXXX9dFF12kfv36adiwYbrnnnsyFRIAwHJ9MtXwtm3bdPz4cf32t7/VmDFj9Oabb2ru3Lk6cuSIfvnLX0qS2traNGXKFE2ePFmrVq3SG2+8oZ/85CcqLS3VNddck6nQAACWynMcx8nWyZYvX66VK1fq3XfflSStXLlSP//5z9Xc3Ky+fftKkhYtWqQnnnhC27Zt67SNeDyueDye2G5tbdXw4cN1of5BfVSQ+Q8BAEiro+rQS3pGhw4dUklJSUptZexOrTOtra0qKytLbDc0NOg73/lOIqFJ0tSpU3X33Xfr448/1sCBA5PaqK+v15IlS5J+/5KeyUzQAICsOHDgQO9Jajt27NADDzyQ6HqUpObmZo0aNcpz3ODBgxP7OktqdXV1qq2tTWwfOnRII0aMUFNTU8oXI5va2to0bNgw7dmzR8XFxVGH0229NW6p98ZO3NlF3Nn39x43901PT4VOaosWLdLdd9/te8xf//pXnX766Ynt999/X9OmTdPll1+uuXPnho/SJRaLKRaLJf2+pKSk132RklRcXEzcWdZbYyfu7CLu7MvPT/3ZxdBJbeHChZozZ47vMaNHj078/MEHH2jSpEk6//zz9W//9m+e4yorK9XS0uL53d+3Kysrw4YGADjBhU5q5eXlKi8v79ax77//viZNmqQJEyZo9erVSVm4urpaP//5z9XR0aGCgi8e8nj++ec1duzYTrseAQDwk7Fxau+//74mTpyo4cOH65e//KU+/PBDNTc3q7m5OXHMD3/4Q/Xt21dXX3213nrrLT366KP613/9V0/NLEgsFtNtt93WaZdkLiPu7OutsRN3dhF39qUz9ow90r9mzRpdddVVne5zn/L1119XTU2N/vKXv2jQoEG6/vrrdcstt2QiJACA5bI6Tg0AgExi7kcAgDVIagAAa5DUAADWIKkBAKzRa5Nab17a5s4779T555+voqIilZaWdnpMXl5e0mvdunXZDdTQnbibmpp0ySWXqKioSBUVFfrpT3+qo0ePZjfQbhg5cmTS9V22bFnUYSVZsWKFRo4cqX79+qmqqkqvvPJK1CEFuv3225OurXuGoVzx4osv6tJLL9XQoUOVl5enJ554wrPfcRwtXrxYQ4YMUWFhoSZPnqx33nknmmBdguKeM2dO0vWfNm1aNMG61NfX61vf+pYGDBigiooKzZw5U9u3b/cc8/nnn6umpkannHKK+vfvr1mzZiVN0BGk1yY199I2b731lu677z6tWrVKP/vZzxLH/H1pmxEjRmjr1q1avny5br/99qSZTbKtvb1dl19+uebNm+d73OrVq7Vv377Ea+bMmdkJsAtBcR87dkyXXHKJ2tvbtXnzZq1du1Zr1qzR4sWLsxxp9yxdutRzfa+//vqoQ/J49NFHVVtbq9tuu03/8z//o3Hjxmnq1Knav39/1KEF+sY3vuG5ti+99FLUISU5cuSIxo0bpxUrVnS6/5577tGvf/1rrVq1Si+//LJOPvlkTZ06VZ9//nmWI/UKiluSpk2b5rn+jzzySBYj7NymTZtUU1OjLVu26Pnnn1dHR4emTJmiI0eOJI5ZsGCBnnrqKT322GPatGmTPvjgA1122WXhTuRY5J577nFGjRqV2H7wwQedgQMHOvF4PPG7W265xRk7dmwU4SVZvXq1U1JS0uk+Sc7jjz+e1Xi6q6u4n3nmGSc/P99pbm5O/G7lypVOcXGx5zvIBSNGjHDuu+++qMPwdd555zk1NTWJ7WPHjjlDhw516uvrI4wq2G233eaMGzcu6jBCMf+9HT9+3KmsrHSWL1+e+N2hQ4ecWCzmPPLIIxFE2LnO/jsxe/ZsZ8aMGZHEE8b+/fsdSc6mTZscx/ni+hYUFDiPPfZY4pi//vWvjiSnoaGh2+322ju1znR3aZvt27fr448/jiLEUGpqajRo0CCdd955euihhzyD1nNRQ0ODzjrrrMRKC9IX17utrU1vvfVWhJF1btmyZTrllFM0fvx4LV++PKe6Sdvb27V161ZNnjw58bv8/HxNnjxZDQ0NEUbWPe+8846GDh2q0aNH60c/+pGampqiDimUXbt2qbm52XP9S0pKVFVV1Suu/8aNG1VRUaGxY8dq3rx5OnDgQNQhJWltbZWkxH+zt27dqo6ODs81P/300zV8+PBQ1zyr66llUrqWtskVS5cu1Xe/+10VFRXpueee03XXXafDhw/rhhtuiDq0LjU3N3sSmuS93rnkhhtu0LnnnquysjJt3rxZdXV12rdvn+69996oQ5MkffTRRzp27Fin17OrBXRzRVVVldasWaOxY8dq3759WrJkiS666CK9+eabGjBgQNThdcvf/147u/659rdsmjZtmi677DKNGjVKO3fu1M9+9jNNnz5dDQ0NOumkk6IOT5J0/Phx3Xjjjbrgggt05plnSlJisWizXh/2mufcndqiRYs6fUjC/TL/UadzaZtsxu3n1ltv1QUXXKDx48frlltu0c0336zly5fnfNxRCvNZamtrNXHiRJ199tm69tpr9atf/UoPPPCAZ1V19Mz06dN1+eWX6+yzz9bUqVP1zDNfrGj8xz/+MerQTgjf//739Y//+I8666yzNHPmTD399NP6y1/+oo0bN0YdWkJNTY3efPPNjDz8lnN3ar11aZuwcYdVVVWlO+64Q/F4PK0TlqYz7srKyqSn87K5lFAqn6WqqkpHjx7V7t27NXbs2AxEF86gQYN00kkndfr329uWZSotLdVpp52mHTt2RB1Kt/39Gre0tGjIkCGJ37e0tOicc86JKKqeGT16tAYNGqQdO3bo4osvjjoczZ8/X08//bRefPFFffWrX038vrKyUu3t7Tp06JDnbi3s33zOJbXeurRNmLh7orGxUQMHDkz7DNzpjLu6ulp33nmn9u/fr4qKCklfXO/i4mKdccYZaTmHn1Q+S2Njo/Lz8xNxR61v376aMGGC1q9fn3jq9fjx41q/fr3mz58fbXAhHT58WDt37tSVV14ZdSjdNmrUKFVWVmr9+vWJJNbW1qaXX3458KnlXLN3714dOHDAk5yj4DiOrr/+ej3++OPauHFjUmlowoQJKigo0Pr16zVr1ixJ0vbt29XU1KTq6upQJ+qV9u7d64wZM8a5+OKLnb179zr79u1LvP7u0KFDzuDBg50rr7zSefPNN51169Y5RUVFzm9/+9sII3ec9957z3nttdecJUuWOP3793dee+0157XXXnM++eQTx3Ec5z/+4z+cf//3f3feeOMN55133nEefPBBp6ioyFm8eHFOx3306FHnzDPPdKZMmeI0NjY6zz77rFNeXu7U1dVFGrdp8+bNzn333ec0NjY6O3fudP7whz845eXlzo9//OOoQ/NYt26dE4vFnDVr1jhvv/22c8011zilpaWep0tz0cKFC52NGzc6u3btcv77v//bmTx5sjNo0CBn//79UYfm8cknnyT+hiU59957r/Paa6857733nuM4jrNs2TKntLTUefLJJ53XX3/dmTFjhjNq1Cjns88+y9m4P/nkE+emm25yGhoanF27djkvvPCCc+655zpf+9rXnM8//zzSuOfNm+eUlJQ4Gzdu9Pz3+tNPP00cc+211zrDhw93/vznPzuvvvqqU11d7VRXV4c6T69NaqtXr3Ykdfpy+9///V/nwgsvdGKxmPOVr3zFWbZsWUQRf2n27Nmdxr1hwwbHcRznP//zP51zzjnH6d+/v3PyySc748aNc1atWuUcO3Ysp+N2HMfZvXu3M336dKewsNAZNGiQs3DhQqejoyO6oDuxdetWp6qqyikpKXH69evnfP3rX3fuuuuuyP/Rd+aBBx5whg8f7vTt29c577zznC1btkQdUqArrrjCGTJkiNO3b1/nK1/5inPFFVc4O3bsiDqsJBs2bOj073n27NmO43zxWP+tt97qDB482InFYs7FF1/sbN++PdqgHf+4P/30U2fKlClOeXm5U1BQ4IwYMcKZO3duTvyPUFf/vV69enXimM8++8y57rrrnIEDBzpFRUXO9773Pc+NSnew9AwAwBo59/QjAAA9RVIDAFiDpAYAsAZJDQBgDZIaAMAaJDUAgDVIagAAa5DUAADWIKkBAKxBUgMAWIOkBgCwxv8DEHJB30KqDHQAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "xbins, ybins, surf_density = vis.density_projection(s, bins=np.linspace(-20, 20, 100))" ] }, { "cell_type": "markdown", "id": "611b13e0", "metadata": {}, "source": [ "## Misc Analyses\n", "\n", "Other analyses are outlined here." ] }, { "cell_type": "code", "execution_count": 16, "id": "eb2bd0ba", "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$[15306613] \\; \\mathrm{\\frac{M_{\\odot}}{kpc^{3}}}$" ], "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# To find the particle density near a given point,\n", "s.density_points([0,0,0]*u.kpc)" ] }, { "cell_type": "code", "execution_count": 17, "id": "1441075f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([2.56199069e+12, 2.49863973e+12, 2.52594940e+12]),\n", " array([[-0.5373118 , 0.72543306, -0.43016614],\n", " [ 0.77115326, 0.62910866, 0.09769822],\n", " [-0.34149476, 0.27922962, 0.89744758]]))" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# To calculate the principal axes (and their directions) of the \n", "# selected particle distribution by diagonalizing the inertia tensor,\n", "s.principal_axes()" ] } ], "metadata": { "kernelspec": { "display_name": "py311", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.4" } }, "nbformat": 4, "nbformat_minor": 5 }