{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example notebook for using Cantor data and off-centre FOF haloes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In a nutshell: Cantor is a re-processing of the Subfind subhalo catalogue, to more accurately assign particles to subhaloes. It is based on the same FOF groups as the Subfind catalogue." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using the Cantor data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Cantor data are stored in the `.../highlev/` director of each (hydro) simulation, in a subdirectory called `Cantor`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`bahe@watering:/net/quasar/data3/Hydrangea/CE-0/HYDRO/highlev/Cantor>ls\n", "Cantor_000.hdf5 Cantor_008.hdf5 Cantor_016.hdf5 Cantor_024.hdf5\n", "Cantor_001.hdf5 Cantor_009.hdf5 Cantor_017.hdf5 Cantor_025.hdf5\n", "Cantor_002.hdf5 Cantor_010.hdf5 Cantor_018.hdf5 Cantor_026.hdf5\n", "Cantor_003.hdf5 Cantor_011.hdf5 Cantor_019.hdf5 Cantor_027.hdf5\n", "Cantor_004.hdf5 Cantor_012.hdf5 Cantor_020.hdf5 Cantor_028.hdf5\n", "Cantor_005.hdf5 Cantor_013.hdf5 Cantor_021.hdf5 Cantor_029.hdf5\n", "Cantor_006.hdf5 Cantor_014.hdf5 Cantor_022.hdf5 Cantor_029_IDs.hdf5\n", "Cantor_007.hdf5 Cantor_015.hdf5 Cantor_023.hdf5 GalaxyTables.hdf5`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are three kinds of HDF5 files: (i) the actual catalogues with properties of subhaloes (`Cantor_0xx.hdf5`); (ii) a separate file for each catalogue the contains the IDs of each particle within subhaloes (`Cantor_0xx_IDs.hdf5`, currently only copied to Leiden for snapshot 29 because these are quite big); (iii) a file `GalaxyTables.hdf5` that connects each galaxy ID to its Cantor subhalo index in each snapshot." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cantor files in more detail" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Galaxy connection tables (GalaxyTables.hdf5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the simplest file, and there is only one per simulation. It gives the Cantor subhalo index for every galaxy in each snapshot, analogous to the `SHI` data set in `FullGalaxyTables.hdf5`. In addition, there is a data set pointing to the Cantor definition, which is not guaranteed to be the same as the central assigned by Subfind (Cantor uses a regularization procedure to avoid central/satellite swaps):\n", "\n", "`bahe@watering:/net/quasar/data3/Hydrangea/CE-0/HYDRO/highlev/Cantor>h5ls GalaxyTables.hdf5 \n", "CentralGalaxy Dataset {389777, 30}\n", "Header Group\n", "SubhaloIndex Dataset {389777, 30}`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Cantor catalogue files (Cantor_0xx.hdf5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each Cantor catalogue file contains two base groups (and a Header group with mostly technical info): `FOF` and `Subhalo`: " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`bahe@watering:/net/quasar/data3/Hydrangea/CE-0/HYDRO/highlev/Cantor>h5ls Cantor_029.hdf5\n", "FOF Group\n", "Header Group\n", "Subhalo Group`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`FOF` does not hold much info, because the FOF properties are the same as in the Subfind catalogues. So the only info here is on the (Cantor) subhaloes within each FOF:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`bahe@watering:/net/quasar/data3/Hydrangea/CE-0/HYDRO/highlev/Cantor>h5ls Cantor_029.hdf5/FOF\n", "CenSubhalo Dataset {104067} --> pointer to central subhalo in each FOF (see below)\n", "FirstSubhalo Dataset {104068} --> pointer to first subhalo in each FOF (see below)\n", "Length Dataset {104067} --> number of particles in ID list for each FOF group\n", "NumOfSubhaloes Dataset {104067} --> number of subhaloes in each FOF group\n", "Offset Dataset {104068} --> offset of particles in ID list for each FOF group` " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that `CenSubhalo` and `FirstSubhalo` are often identical; they only differ where a FOF has not a single subhalo (`NumOfSubhaloes == 0`). In those cases, `CenSubhalo = -1`, whereas `FirstSubhalo` points to the central subhalo in the *next* FOF that has any (this sounds silly at first, but is useful in certain situations because it makes the list monotonic)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`Subhalo` is (much) more interesting. It contains (a currently limited set of) properties of individual Cantor subhaloes. For most of these, it should be relatively clear from the name what they mean, there is also a more detailed description in the `Comment` attribute of each." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`bahe@watering:/net/quasar/data3/Hydrangea/CE-0/HYDRO/highlev/Cantor>h5ls Cantor_029.hdf5/Subhalo\n", "AngularMomentum_DM Dataset {121058, 3}\n", "AngularMomentum_Gas Dataset {121058, 3}\n", "AngularMomentum_Stars Dataset {121058, 3}\n", "CentreOfMass Dataset {121058, 3}\n", "CentreOfPotential Dataset {121058, 3}\n", "Extra Group\n", "FOF_Index Dataset {121058}\n", "Galaxy Dataset {121058}\n", "IndexBySubfindID Dataset {114234}\n", "Length Dataset {121058}\n", "LengthType Dataset {121058, 6}\n", "Mass Dataset {121058}\n", "MassType Dataset {121058, 6}\n", "MaxRadius Dataset {121058}\n", "MaxRadiusType Dataset {121058, 6}\n", "Offset Dataset {121058}\n", "OffsetType Dataset {121058, 7}\n", "Position Dataset {121058, 3}\n", "RadiusOfVmax Dataset {121058}\n", "SubfindIndex Dataset {121058}\n", "Velocity Dataset {121058, 3}\n", "VelocityDispersion_DM Dataset {121058}\n", "VelocityDispersion_Stars Dataset {121058}\n", "Vmax Dataset {121058}\n", "ZMF_Velocity Dataset {121058, 3}`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Example for more detailed description (in this case for the `LengthType` data set):\n", "\n", "`bahe@watering:/net/quasar/data3/Hydrangea/CE-0/HYDRO/highlev/Cantor>h5dump -A -d Subhalo/LengthType Cantor_029.hdf5\n", "HDF5 \"Cantor_029.hdf5\" {\n", "DATASET \"Subhalo/LengthType\" {\n", " DATATYPE H5T_STD_I32LE\n", " DATASPACE SIMPLE { ( 121058, 6 ) / ( 121058, 6 ) }\n", " ATTRIBUTE \"Comment\" {\n", " DATATYPE H5T_STRING {\n", " STRSIZE 164;\n", " STRPAD H5T_STR_NULLPAD;\n", " CSET H5T_CSET_ASCII;\n", " CTYPE H5T_C_S1;\n", " }\n", " DATASPACE SCALAR\n", " DATA {\n", " (0): \"Number of particles in ID list that belong to subhalo i (first index) and have type j (second index). These particles are stored at indices [offset_i]:[offset_i+1].\"\n", " }\n", " }\n", "}\n", "}`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Maybe worth highlighting particularly is the `IndexBySubfindID` data set: that one contains the Cantor subhalo index for a given *Subfind* subhalo index in the same snapshot (-1 if that subfind subhalo was not found by Cantor):" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/software/rhel7/lib64/python3.6/site-packages/dask/config.py:168: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.\n", " data = yaml.load(f.read()) or {}\n" ] } ], "source": [ "import hydrangea as hy\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/net/quasar/data3/Hydrangea/10r200/CE-0/HYDRO/highlev/Cantor/Cantor_029.hdf5\n" ] } ], "source": [ "sim = hy.Simulation(0) # Set up simulation object for CE-0, as an example\n", "cantor_cat = f'{sim.high_level_dir}Cantor/Cantor_029.hdf5'\n", "print(cantor_cat)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0 1 2 3 4 5 6 7 8 10 9 12 13 11 17 14 18 32 19 16]\n" ] } ], "source": [ "cantor_for_subfind = hy.hdf5.read_data(cantor_cat, 'Subhalo/IndexBySubfindID')\n", "print(cantor_for_subfind[:20])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So here, for example, Subfind subhalo 9 corresponds to Cantor subhalo 10 (and, by coincidence, vice versa). We can directly compare e.g. the coordinates and stellar mass of this object in both catalogues:" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Coordinates in Subfind are [1608.3173 1602.9695 1607.5225]\n", "Coordinates in Cantor are [1608.31688088 1602.96922824 1607.52246131]\n", "\n", "Stellar mass in Subfind is 1.1915 x 10^11 M_Sun.\n", "Stellar mass in Cantor is 1.2450 x 10^11 M_Sun.\n" ] } ], "source": [ "subfind_file = sim.get_subfind_file(29)\n", "sub = hy.SplitFile(subfind_file, 'Subhalo', read_index=9, verbose=0)\n", "print(\"Coordinates in Subfind are\", sub.CentreOfPotential)\n", "centre = hy.hdf5.read_data(cantor_cat, 'Subhalo/CentreOfPotential', read_index=10)\n", "print(\"Coordinates in Cantor are\", centre)\n", "print(\"\")\n", "print(f\"Stellar mass in Subfind is {sub.MassType[4]/1e11:.4f} x 10^11 M_Sun.\")\n", "print(f\"Stellar mass in Cantor is \"\n", " f\"{hy.hdf5.read_data(cantor_cat, 'Subhalo/MassType', read_index=10)[4]/10:.4f} x 10^11 M_Sun.\")\n", " \n", "## Careful: masses in the Cantor catalogue are in 10^10 M_Sun, but those extracted from the Subfind catalogue\n", "## with the SplitFile class are in M_Sun..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Result: the coordinates are almost exactly identical (the tiny difference comes from identifying a slightly different particle as the potential minimum), but the stellar masses are not quite the same. Cantor finds a bit of a higher stellar mass, because this one is a satellite (the central cluster, `FOF == 0` is bound to have at least 10 subhaloes!), and Cantor tends to attach more stars to a subhalo than Subfind.\n", "\n", "Also note that the Cantor coordinates are stored at a higher precision, which is important since the Subfind (float) precision only allows ~kpc accuracy." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Also worth pointing out: the `Galaxy` data set points to the galaxy ID of each Cantor subhalo. This is the same galaxy ID as in the `FullGalaxyTables.hdf5` catalogues (since Cantor is built on top of these galaxy tracing results). Let's test this explicitly:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Galaxy ID of Cantor subhalo 10 is 767.\n", "Subfind subhalo of galaxy 767 in snapshot 29 is 9.\n" ] } ], "source": [ "galaxy = hy.hdf5.read_data(cantor_cat, 'Subhalo/Galaxy', read_index=10)\n", "print(f\"Galaxy ID of Cantor subhalo 10 is {galaxy}.\")\n", " \n", "# Find the Subfind subhalo of this one from the FullGalaxyTables file\n", "subfind_from_galaxy = hy.hdf5.read_data(sim.fgt_loc, 'SHI', read_index=galaxy)[29] # Look up subhalo in snap 29\n", "print(f\"Subfind subhalo of galaxy {galaxy} in snapshot 29 is {subfind_from_galaxy}.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, the `Extra` group (within `Subhalo`) contains some more detailed subhalo properties that are only calculated for reasonably massive subhaloes. There are two datasets within `Extra` to find those subhaloes: (i) `ExtraIDs` gives the index *in the reduced \"extra properties\" catalogue* for each (Cantor) subhalo; for all subhaloes that are not massive enough, this is `-1`. (ii) `SubhaloIndex` gives the reverse: the index *in the full (Cantor) subhalo catalogue* for each entry in the extra catalogue.\n", "\n", "As an example, again with our subhalo 10 from above:" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Subhalo 10 has extra-ID 10.\n" ] } ], "source": [ "extraID = hy.hdf5.read_data(cantor_cat, 'Subhalo/Extra/ExtraIDs', read_index=10)\n", "print(f\"Subhalo 10 has extra-ID {extraID}.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So in this case, because the top 10 subhaloes are all massive enough to be in the extra catalogue, the two indices are the same." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Linking to particles (Cantor_0xx_IDs.hdf5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To find all particles that Cantor identifies as belonging to a particular subhalo, you can use the `Cantor_0xx_IDs.hdf5` files (e.g. `Cantor_029_IDs.hdf5` for z = 0). This contains one single data set `IDs`:\n", "\n", "`bahe@watering:/net/quasar/data3/Hydrangea/CE-0/HYDRO/highlev/Cantor>h5ls Cantor_029_IDs.hdf5\n", "IDs Dataset {109134215}`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This works very similar to Subfind. `IDs` contains the particle IDs of each Cantor subhalo, in order of both FOF groups and Cantor subhaloes (this works because subhaloes are labelled sequentially within each FOF group). In addition, within each subhalo, IDs are sorted by particle type, and within particle type by radial distance from the subhalo centre.\n", "\n", "The information about where in this long list the IDs for a particular subhalo are sitting is specified in the `Offset...` and `Length...` data sets in the main Cantor catalogue file (`Cantor_0xx.hdf5`, as described above). There are three of these; see below for examples:\n", "\n", "(i) `Offset`: index (in the `IDs` data set) of the *first* particle belonging to each subhalo;\n", " `Length`: total number of particles belonging to this subhalo. \n", " \n", "(ii) `OffsetType` (and `LengthType`): the same, but with a second dimension for each particle type, so it specifies the block of e.g. gas particles of a given subhalo.\n", "\n", "(iii) `Extra/OffsetTypeApertures`: the same, but with a *third* dimension for different radial apertures. So this can be used to directly find particles of a given type that lie within a certain radius from the subhalo centre (3, 10, 30, and 100 pkpc)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Example 1: whole subhalo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's retrieve a list with the IDs of all particles belonging to our Cantor subhalo 10 from above:" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Retrieved 127373 IDs.\n", "First 10 IDs are: [ 90454593 166190593 166177729 335199581 204404587 240929347 136430095\n", " 310787423 282912803 324814189]\n" ] } ], "source": [ "offset = hy.hdf5.read_data(cantor_cat, \"Subhalo/Offset\", read_index=10)\n", "length = hy.hdf5.read_data(cantor_cat, \"Subhalo/Length\", read_index=10)\n", "\n", "# Find the first particle *after* the block belonging to our subhalo\n", "# Note that we cannot (reliably) use the offset of the next subhalo (11 here), because \n", "# at the end of the particle block for each FOF group there are \"loose\" particles that\n", "# do not belong to any subhalo.\n", "end = offset + length\n", "\n", "cantor_ids = f'{sim.high_level_dir}Cantor/Cantor_029_IDs.hdf5'\n", "ids = hy.hdf5.read_data(cantor_ids, 'IDs', read_range=[offset, end])\n", "print(f\"Retrieved {ids.shape[0]} IDs.\")\n", "print(\"First 10 IDs are:\", ids[:10])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Example 2: only star particles" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In example 1, we got *all* particles. Let's do this a bit more finely now, only retrieving stars:" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Retrieved 84120 IDs.\n", "First 10 IDs are: [348327033 151511633 50898603 24679775 209739697 348974305 271038371\n", " 275338903 44062981 255913957]\n" ] } ], "source": [ "offset = hy.hdf5.read_data(cantor_cat, \"Subhalo/OffsetType\", read_index=10)[4]\n", "end = hy.hdf5.read_data(cantor_cat, \"Subhalo/OffsetType\", read_index=10)[5]\n", "\n", "# Note: in this case we *can* directly look up the offset of the next particle type as `end`,\n", "# because these are subdivisions of individual subhaloes; there is an extra \"coda\"\n", "# entry in index 6 of the 2nd dimension so this also works for BHs (type 5).\n", "\n", "ids_stars = hy.hdf5.read_data(cantor_ids, 'IDs', read_range=[offset, end])\n", "print(f\"Retrieved {ids_stars.shape[0]} IDs.\")\n", "print(\"First 10 IDs are:\", ids_stars[:10])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This list is different from above, because the full ID list starts with the gas particles (type 0). We can verify that these are indeed in the full list though:" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True\n", "True\n", "True\n", "True\n", "True\n", "True\n", "True\n", "True\n", "True\n", "True\n" ] } ], "source": [ "for id in ids_stars[:10]:\n", " print(id in ids)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To test that the result is sensible, plot a histogram of star particle radii:" ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Furthest star in subhalo 10 is at r = 0.767 Mpc\n", "Checking 8 cells...\n", "Region setup took 0.010 sec.\n", "Selection region contains 8 cells, 11 segments, 1779637 particles, 4 files\n", "Reading 'ParticleIDs' took 0.148 sec.\n", "Located 84120 out of 84120 particles in snapshot region.\n", "Reading 'Coordinates' took 0.233 sec.\n", "Prepared reading from 'IDS'...\n", "Loaded file offsets in 0.159 sec.\n", "46 \n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAArb0lEQVR4nO3deZxU1ZXA8d+RHUFZbJEADhpwIW5gq6hN44aicQABE9BRNERM1EnUMeOSTPjERINOosYYzaComDEKrqCiDSFoBAVphVZxo8Ug9LC0sigqIHDmj3NLCuimq7ur6r2qOt/Ppz716tarqlMFferVuffdK6qKc865wrBH1AE455zLHk/6zjlXQDzpO+dcAfGk75xzBcSTvnPOFZCmUQewO/vss49279496jCccy6nvP7665+oalFN98U66Xfv3p3y8vKow3DOuZwiIktru8/LO845V0A86TvnXAHxpO+ccwUkpaQvIleJyCIReVtEHhGRliJygIjME5FKEZkkIs3Dvi3C7cpwf/ek57k+tL8vImdk6D0555yrRZ1JX0S6AD8BilX1MKAJMAK4BbhdVXsAa4HR4SGjgbWh/fawHyLSKzzuO8BA4G4RaZLet+Occ253Ui3vNAVaiUhToDWwAjgFeDzcPxEYErYHh9uE+08VEQntj6rqJlX9CKgEjm30O3DOOZeyOpO+qlYBvwM+xpL9euB1YJ2qbgm7LQe6hO0uwLLw2C1h/47J7TU85hsiMkZEykWkvLq6uiHvyTnnXC1SKe+0x47SDwC+BeyJlWcyQlXHq2qxqhYXFdV4boFzWfPZZ/Dgg7B5c9SROJceqZR3TgM+UtVqVf0aeBI4EWgXyj0AXYGqsF0FdAMI9+8NfJrcXsNjXB7Ytg3GjoW77oo6kvTYvBmGDYOLL4Z77406GufSI5Wk/zHQV0Rah9r8qcA7wCxgeNhnFDAlbE8Ntwn3/11tpZapwIgwuucAoCfwWnrehovali1w0UVw441w/fXw1VdRR9Q4qnDJJfC3v8F++8Ftt8HWrVFH5VzjpVLTn4d1yL4BvBUeMx64FrhaRCqxmv2E8JAJQMfQfjVwXXieRcBk7AvjBeByVfU/ozyweTOMHAl/+YsdGW/YAM8+G3VUjTN2LDz0EPzqV/bLZckSeOqpqKNyrvEkzsslFhcXq8+9E28bN8Lw4fDcc3Y0/JOfQLdu0LcvPPlk1NE1zH332VH+6NFW1tm2DQ4+GDp2hLlzQSTqCJ3bPRF5XVWLa7rPz8h1DfbFF3D22TBtGvz5z3DVVdCkCYwYYV8C69ZFHWH9vfAC/OhHcMYZcM89luCbNIGrr4bXXoPZs6OO0LnG8aTvGmT9ekuMs2bBxIlw6aXb7xs50ko+uXakv2ABnHsuHHEEPPYYNGu2/b6LLrIj/d/9LrLwnEsLT/qu3tasgdNOg3nz4NFH4YILdry/uBh69IC//jWa+Bpi6VI46yzo0MH6I9q23fH+1q3h8sth6lR4771oYnQuHTzpu3pZtQpOOgneess6Ns89d9d9ROC88+xXwIoVWQ+x3tauhTPPtBFH06bBt75V836XXw4tW1rfhXO5ypO+S9ny5VBaCh9+aEfDZ59d+74jR1oH6OTJ2YuvITZtgnPOsff09NPwne/Uvu+++8KoUTaqZ9WqrIXoXFp50ncp+egjS/grVkBZmZV3dueQQ6B373iXeLZts1r9Sy/ZWbcnnVT3Y66+2vor8uUENFd4POm7On3yCfTrZ6NxZs6EkpLUHnfeeTbipbIyo+E12PXXW5/EuHH2yyQVBx0EgwfD3Xfb6CXnco0nfVenGTOgqspq+Mcck/rjRoyw+v4jj2Qutob6n/+BW2+FH/8Y/vM/6/fYa66xzuwHHshMbM5lkid9V6eKChu+ePzx9Xtc165WEvrrX21ag7hYt84S/YABcOed9T/Z6oQT7OQzn5rB5SJP+q5OFRXQqxc0b17/x44caUMcKyrSH1dD/elPNnvmLbdA06Z1778zEfjZz6yfI9fORXDOk76rU0UFHHlkwx47fLgl1rh06H7xBdxxh43J79274c8zeLCdi/Df/x2vXzHO1cWTvtut1attxE5Dk37HjjBwoNX1t21Lb2wNce+91jF9ww2Ne57E1Azz58PLL6cnNueywZO+261EWaahSR9sFM/y5dHPW7Npk02j0L8/nHhi459v1CjYZx+fmsHlFk/6brfSkfQHDbJpDKIexfPQQzYKqbFH+QmJqRmeeQbefTc9z+lcpnnSd7tVUQFdutgRbUPtuafVwCdPjm7ZwS1bbDx+cbGN2kkXn5rB5RpP+m63GtOJm+y882xs+4wZjX+uhpg82RZC+fnP0zsfflGRndX70EOwcmX6nte5TEllYfSDRWRh0uUzEblSRDqIyAwRWRyu24f9RUTuFJFKEXlTRPokPdeosP9iERlV+6u6ONi0ycoW6Uj6p59uM1hGMYpn2za4+WabV2fQoPQ//9VXw9df+9QMLjekslzi+6p6lKoeBRwNfAk8hS2DOFNVewIzw22AM7H1b3sCY4B7AESkAzAWOA44Fhib+KJw8fTuu1YWSUfSb97chm9OmZL96QueeQYWLbJpF/bIwG/bnj1hyBCbmmHDhvQ/v3PpVN8/gVOBD1V1KTAYmBjaJwJDwvZg4CE1c4F2ItIZOAOYoaprVHUtMAMY2Ng34DInHZ24yc47zxL+M8+k5/lSoQo33QQHHgjf/37mXueaa2yKZp+awcVdfZP+CCAxBqOTqiZmS18JdArbXYBlSY9ZHtpqa9+BiIwRkXIRKa+urq5neC6dKiqgVSs7kk2Hfv2sUzibJZ6ZM20s/bXXNuzs21SdcIJ1Ev/v/2buNZxLh5STvog0BwYBj+18n9rq6mk5L1FVx6tqsaoWFxUVpeMpXQNVVMBhh9mJSOmwxx42LcPzz8Onn6bnOety0022KMqoLPQgDRgAb7zhs2+6eKvPkf6ZwBuqmlg+YlUo2xCuV4f2KqBb0uO6hrba2l0MqaZv5E6ykSOtn+CJJ9L7vDV55RV48UUrvbRokfnXKymx9zZvXuZfy7mGqk/SH8n20g7AVCBx/DQKmJLUfmEYxdMXWB/KQGXA6SLSPnTgnh7aXAxVVdnReLqTfu/ecPDB2TlR6+abbRqIMWMy/1pgJR6R6M88dm53Ukr6IrInMABInlNwHDBARBYDp4XbANOAJUAlcC9wGYCqrgF+DcwPlxtDm4uhRCfuUUel93kT6+e+9JJNzZApCxfCc8/BlVfayWHZ0K4dHH64z8Xj4i2lpK+qX6hqR1Vdn9T2qaqeqqo9VfW0RAIPo3YuV9Vvq+rhqlqe9Jj7VbVHuPg4hxhLJP0jjkj/c48caeWjSZPS/9wJv/0ttG1rZ8xmU79+8OqrVuZxLo78jFxXo4oKOOAA2Guv9D93z5420iVTo3g++AAee8wSfvssnwlSUmIduXFaP8C5ZJ70XY0y0Ymb7LzzbKTL/Pnpf+5x46zj9sor0//cdUmsH+wlHhdXnvTdLr78EhYvzmzSv/BCW05x2LD0zlmzdCn85S9wySXQqVPd+6db167Qvbt35rr48qTvdvH22zZfTSaTfseOMHWqjRA65xzYuDE9z5uY2/6aa9LzfA1RUmJJ31fUcnHkSd/tIt3TL9Smd287Kp87F374w8YnyY8+gvvus18R+++fnhgboqQEVq2CysroYnCuNp703S4qKmzkS/fumX+toUPhN7+Bhx+2WnxDvfmmJdsWLWxitSj162fXXuJxceRJ3+2iosKGamZiRsqa3HCDdezecAM89VT9Hz9rliVaEetA7dEj/THWxyGH2DTSnvRdHHnSdzvYti3zI3d2JmJlmWOPhX/7NzuxKlWTJtnC61272vj4ww/PWJgp22MPW4PXR/C4OPKk73bwz3/C55+n/0zcurRqBU8/bUfIgwalNqLnjjtgxAj7snj5ZejWrc6HZE2/fjYCatWquvd1Lps86bsdZKsTtyadO9siK598svsRPdu22eicq66yPoEZM+zLIk4S4/XnzIk2Dud25knf7aCiwsoThx0Wzev36bN9RM8ll+w6omfzZrjgAvj97+2M28mTbWHyuDn6aIvLSzwubjzpux1UVNg0Ca1bRxfDsGHw61/bgiTJI3o++wzOOsumb7j5ZvjjH9M313+6NW8Oxx3nnbkufjK4lpDLRRUVcMwxUUcBP/85vPOOjeg59FBLoGeeaWvdPvhgdhZFaaySEvvS2rAB2rSJOhrnjB/pu2989pmd4BRFPX9nIjBhgn0BnX8+9O1rJzs980xuJHywpL91q5WqnIsLT/ruG2++addxSPpgI3qmTLGZMjdutFWwBg6MOqrUnXCC9Y94icfFiZd33DeiHLlTm86dLS6R+I3Qqctee9lJbt6Z6+Ik1ZWz2onI4yLynoi8KyLHi0gHEZkhIovDdfuwr4jInSJSKSJvikifpOcZFfZfLCI58iO9cFRUWGLt0iXqSHbUsWPuJfyEkhIr73z9ddSROGdSLe/8AXhBVQ8BjgTeBa4DZqpqT2BmuA22gHrPcBkD3AMgIh2AscBxwLHA2MQXhYuHhQvtKF8k6kjyR79+NlV1fc4ydi6T6kz6IrI3UApMAFDVzaq6DhgMTAy7TQSGhO3BwENh2cS5QDsR6QycAcxQ1TWquhaYAeRQhTa/bd1qUyrHqbSTD3xRFRc3qRzpHwBUAw+IyAIRuS8slN5JVVeEfVYCiSUrugDLkh6/PLTV1r4DERkjIuUiUl5dXV2/d+MabPFi+Oqr7E+/kO++9S048EDvzHXxkUrSbwr0Ae5R1d7AF2wv5QC2GDqQliUjVHW8qharanFRUVE6ntKlII6duPnCF1VxcZJK0l8OLFfVeeH249iXwKpQtiFcrw73VwHJU191DW21tbsYqKiApk3tRCiXXiUlUF1tC7Y7F7U6k76qrgSWicjBoelU4B1gKpAYgTMKmBK2pwIXhlE8fYH1oQxUBpwuIu1DB+7poc3FQEWFJfwWLaKOJP/4oiouTlIdp//vwMMi0hxYAlyMfWFMFpHRwFLge2HfacBZQCXwZdgXVV0jIr8G5of9blTVNWl5F67RKirg5JOjjiI/HXywDTudPRtGj446GlfoUkr6qroQKK7hrlNr2FeBy2t5nvuB++sRn8uCTz+Fqiqv52eKiJV4fASPiwOfhsF5J24W9OsHH34IK1bUva9zmeRJ33nSz4LEeH2v67uoedJ3LFwI++0H++4bdST5q3dvm0DOk76Lmid9l/WF0AtR8+Y2PbQnfRc1T/oFbvNmW6zEz8TNvJIS+1X12WdRR+IKmSf9AvfeezYDpB/pZ15JiS3q7ouquCh50i9w3ombPccf74uquOh50i9wFRV2Fu5BB0UdSf5r29bKaD5e30XJk36Bq6iAww6zeXdc5pWUwLx51pfiXBQ86RcwVR+5k239+tkU1gsWRB2JK1Se9AvYypU2+6Mn/ezxRVVc1DzpFzDvxM2+/faDHj28M9dFx5N+AUus23rEEZGGUXB8URUXJU/6BayiAvbfH9r78vRZVVJiM5u+917UkbhC5Em/gL32Ghx9dNRRFJ7SUrv+xz+ijcMVJk/6BWrFCliyZHvHosueHj2stu9J30UhpaQvIv8UkbdEZKGIlIe2DiIyQ0QWh+v2oV1E5E4RqRSRN0WkT9LzjAr7LxaRUbW9nsu8OXPs+sQTo42jEInY0f5LL3ld32VffY70T1bVo1Q1sYLWdcBMVe0JzAy3Ac4EeobLGOAesC8JYCxwHHAsMDbxReGyb84cm+q3d++oIylMpaW2WtlHH0UdiSs0jSnvDAYmhu2JwJCk9ofUzAXaiUhn4AxghqquUdW1wAxgYCNe3zXC7Nlw7LE25a/Lvv797dpLPC7bUk36CkwXkddFZExo66SqicXfVgKdwnYXYFnSY5eHttradyAiY0SkXETKq6urUwzP1ccXX9gZoV7Pj06vXtChgyd9l32pzrhSoqpVIrIvMENEdhhspqoqImmpTqrqeGA8QHFxsVc8M2DePNi61ev5UdpjD5uSwZO+y7aUjvRVtSpcrwaewmryq0LZhnC9OuxeBXRLenjX0FZbu8uyOXOsM/H446OOpLD172+LpVf5X4HLojqTvojsKSJtE9vA6cDbwFQgMQJnFDAlbE8FLgyjePoC60MZqAw4XUTahw7c00Oby7LZs21mzXbtoo6ksPl4fReFVI70OwGzRaQCeA14TlVfAMYBA0RkMXBauA0wDVgCVAL3ApcBqOoa4NfA/HC5MbS5LNq6FV591ev5cXDkkTbHvid9l0111vRVdQmwy5RcqvopcGoN7QpcXstz3Q/cX/8wXbq89RZ8/rnX8+OgaVP78vWk77LJz8gtMImTsvxIPx5KS21heh+o5rLFk36BmT0bunSxidZc9BJ1fZ9f32WLJ/0CM2eOlXZEoo7EARQX25nRXuJx2eJJv4B8/DEsW+alnThp3tyGzr70UtSRuELhSb+A+CRr8VRaamsbrFsXdSSuEHjSLyCzZ0ObNr5SVtz072+zbSa+lJ3LJE/6BWTOHOjb14YKuvg47jho1szr+i47POkXiPXrbYy+1/Pjp1Urm/HU6/ouGzzpF4i5c2HbNq/nx1VpKbz+OmzYEHUkLt950i8Qc+ZAkyZWSnDx078/bNliX87OZZIn/QIxe/b2uV5c/Jxwgk237CUel2me9AvA11/bHPpez4+vtm2hTx/vzHWZ50m/ACxcCF9+6fX8uOvf376cN26MOhKXzzzpFwA/KSs3lJbCpk0wf37Ukbh85km/AMyeDd2720RrLr5KSmxOJK/ru0zypJ/nEmd6ej0//jp0gMMP97q+y6yUk76INBGRBSLybLh9gIjME5FKEZkkIs1De4twuzLc3z3pOa4P7e+LyBlpfzduF0uWwMqVXtrJFaWl8Mor1vnuXCbU50j/p8C7SbdvAW5X1R7AWmB0aB8NrA3tt4f9EJFewAjgO8BA4G4RadK48F1dfNGU3FJaCl98AQsWRB2Jy1cpJX0R6Qp8F7gv3BbgFODxsMtEYEjYHhxuE+4/New/GHhUVTep6kfYGrrHpuE9uN2YPdsWQO/VK+pIXCoSi6p4Xd9lSqpH+ncA/wlsC7c7AutUdUu4vRxIdBN2AZYBhPvXh/2/aa/hMd8QkTEiUi4i5dW+hlyjzZmz/cQfF3+dOsHBB3td32VOnalARM4GVqvq61mIB1Udr6rFqlpcVFSUjZfMW2vW2PqrXs/PLaWltnzi1q1RR+LyUSrHfycCg0Tkn8CjWFnnD0A7EUlM0tsVqArbVUA3gHD/3sCnye01PMZlwCuv2LXX83NL//7bZ0V1Lt3qTPqqer2qdlXV7lhH7N9V9XxgFjA87DYKmBK2p4bbhPv/rqoa2keE0T0HAD2B19L2TtwuZs+2edqPOSbqSFx9JOr6XuJxmdCYSu+1wNUiUonV7CeE9glAx9B+NXAdgKouAiYD7wAvAJerqv+AzaA5c+Doo22+dpc7unWzk+k86btMqNcaSqr6IvBi2F5CDaNvVHUjcG4tj78JuKm+Qbr6S5zOf8UVUUfiGqJ/f5g2zU6uE4k6GpdPfExHnnr9dUv8Xs/PTaWlUF0N770XdSQu33jSz1OzZ9v1CSdEG4drGK/ru0zxpJ+n5syBgw6CffeNOhLXEN/+NnTu7EnfpZ8n/TyUmGTNx+fnLhGr67/0kv17OpcunvTz0Pvvw6efetLPdaWlUFUFH30UdSQun3jSz0OJer534uY2r+u7TPCkn4fmzIF99rGavstdvXrZv6NPvubSyZN+Hpo920o7Pr47t4nAaafZeH2fh8eliyf9PLNqFVRWej0/XwwdCqtXby/ZOddYnvTzzIsv2nW/fpGG4dLkzDOhZUt44omoI3H5wpN+nikrs0VTioujjsSlQ5s2MHAgPPkkbNtW9/7O1cWTfh5RhenTrQ7ctF6zKrk4GzbMhm7Onx91JC4feNLPI++8Y8nhDF9yPq+cfbZNke0lHpcOnvTzSFmZXXvSzy/t2sGpp1rS97NzXWN50s8jZWVw6KE2H7vLL8OGwZIlUFERdSQu13nSzxNffWVnbvpRfn4aPNgWt/cSj2usVBZGbykir4lIhYgsEpFfhfYDRGSeiFSKyCQRaR7aW4TbleH+7knPdX1of19EPD2l0T/+ARs3wumnRx2Jy4SiIpuA7ckno47E5bpUjvQ3Aaeo6pHAUcBAEekL3ALcrqo9gLXA6LD/aGBtaL897IeI9MLW2P0OMBC4W0SapPG9FLTp06FFC0sMLj8NHWqd9b6wimuMVBZGV1XdEG42CxcFTgEeD+0TgSFhe3C4Tbj/VBGR0P6oqm5S1Y+ASmpYbtE1TFmZnZDVunXUkbhMOeccu/YSj2uMlGr6ItJERBYCq4EZwIfAOlXdEnZZDnQJ212AZQDh/vXYwunftNfwmOTXGiMi5SJSXl1dXe83VIiWL4dFi7yen++6dIHjj/ek7xonpaSvqltV9SigK3Z0fkimAlLV8aparKrFRUVFmXqZvDJ9ul170s9/w4bBggU2kse5hqjX6B1VXQfMAo4H2olI4rzPrkBV2K4CugGE+/cGPk1ur+ExrhHKymxpvcMOizoSl2lDh9r1U09FG4fLXamM3ikSkXZhuxUwAHgXS/7Dw26jgClhe2q4Tbj/76qqoX1EGN1zANATeC1N76Ngbd0KM2bYqB2fSjn/HXAA9O7tJR7XcKkc6XcGZonIm8B8YIaqPgtcC1wtIpVYzX5C2H8C0DG0Xw1cB6Cqi4DJwDvAC8DlquqzhDdSeTmsXeulnUIybBi8+qpNueFcfdU5LZeqvgn0rqF9CTWMvlHVjcC5tTzXTcBN9Q/T1Wb6dDvCHzAg6khctgwbBr/4hZV4rrgi6mhcrvEzcnNcWRkcfbQtq+cKwyGH2FKKXuJxDeFJP4etXw9z53pppxANG2ZnYfuoZldfnvRz2MyZ1pHrUy8UnqFDbVGVKVPq3te5ZJ70c1hZGbRtayfsuMJy5JFw4IFe4nH150k/R6la0j/lFFtgwxUWESvxzJwJ69ZFHY3LJZ70c9TixbB0qdfzC9mwYfD11/DMM1FH4nKJJ/0c5atkuWOOga5dvcTj6seTfo4qK4MePayu6wrTHnvYzJtlZbBhQ937Owee9HPSpk0wa5aP2nFW4tm4EZ5/PupIXK7wpJ+D5syBL7/00o6DkhLYd18v8bjUedLPQWVl0LQpnHxy1JG4qDVpAkOGwHPP2RG/c3XxpJ+Dpk+HE0+0MfrODRtmNf3EugrO7Y4n/RyzahUsXOilHbfdSSdBu3a+aLpLjSf9HOOrZLmdNW8OgwbB1Kk2bt+53fGkn2PKyqCoCI46KupIXJwMH27rKrzwQtSRuLjzpJ9Dtm2zI/0BA2yMtnMJAwfCfvvBvfdGHYmLu1SWS+wmIrNE5B0RWSQiPw3tHURkhogsDtftQ7uIyJ0iUikib4pIn6TnGhX2Xywio2p7TVezhQttKl0v7bidNWsGP/iBjeJZvjzqaFycpXK8uAX4D1XtBfQFLheRXtgyiDNVtScwM9wGOBNb/7YnMAa4B+xLAhgLHIetuDU28UXhUpOYesFPynI1GT3afg3ef3/Ukbg4qzPpq+oKVX0jbH+OLYreBRgMTAy7TQSGhO3BwENq5gLtRKQzcAa2vu4aVV0LzAAGpvPN5Lvp021K3f32izoSF0cHHmilvwkTbJ0F52pSr8qwiHTH1sudB3RS1RXhrpVAp7DdBViW9LDloa229p1fY4yIlItIebUvC/SNDRvsTFw/yne7M2YMfPyxj9l3tUs56YtIG+AJ4EpV/Sz5PlVVQNMRkKqOV9ViVS0uKipKx1PmhVmzbDie1/Pd7gwaZNMyjB8fdSQurlJK+iLSDEv4D6tq4hSQVaFsQ7heHdqrgG5JD+8a2mprdykoK4PWrW2uFedq07w5XHyxzbH/f/8XdTQujlIZvSPABOBdVb0t6a6pQGIEzihgSlL7hWEUT19gfSgDlQGni0j70IF7emhzddiwAR55xI7yW7SIOhoXdz/8odX0H3gg6khcHKVypH8icAFwiogsDJezgHHAABFZDJwWbgNMA5YAlcC9wGUAqroG+DUwP1xuDG2uDuPHw5o18LOfRR2JywU9etgymvfdZ6N5nEsmVo6Pp+LiYi0vL486jEht2mSjMg46yOr6zqVi0iQYMcLKgt75X3hE5HVVLa7pPj+vM+Yeeshqs9dfH3UkLpcMGQL77OMdum5XnvRjbMsWuOUWOPpoG3/tXKpatICLLoIpU2DlyqijcXHiST/GHn8cPvzQjvJFoo7G5ZpLLrEDhwcfjDoSFyee9GNKFcaNg0MOscWvnauvgw6C/v1tEjbv0HUJnvRj6vnnoaICrr3WZ9R0DTdmDCxZ4oMA3HaeTmLq5puhWzc477yoI3G5bOhQ6NDBO3Tddp70Y+jll22enZ/9zM6wdK6hWraEUaPgqadg9eq693f5z5N+DN18s62ONXp01JG4fHDJJTZv08SJde/r8p8n/ZhZsMCWvLvySptrx7nGOvRQm7Np/HgbIOAKmyf9mBk3Dtq2hcsuizoSl0/GjIHKSnjxxagjcVHzpB8jH3wAjz0Gl18O7dpFHY3LJ8OH2/8pX0PXedKPkVtvtTMpr7wy6khcvmnVCi68EJ54Aj75JOpoXJQ86cfE8uU2z87o0dCpU937O1dfl1wCmzfb/zNXuDzpx8Tvf29nTV5zTdSRuHx12GFwwgneoVvoPOnHwCef2B/i+edD9+5RR+Py2SWXwPvv27kgrjB50o+BO++EL7+0KRecy6TvfQ/23huuugoWLow6GheFVJZLvF9EVovI20ltHURkhogsDtftQ7uIyJ0iUikib4pIn6THjAr7LxaRUTW9ViH6/HP44x9tUrVevaKOxuW71q1tBM/SpdCnD1x6qZ+pW2hSOdJ/EBi4U9t1wExV7QnMDLcBzgR6hssY4B6wLwlgLHAccCwwNvFFUej+/GdYt84XSXHZc+65sHgx/OQnMGEC9OxpfUqbN0cdmcuGOpO+qv4D2Hkt28FA4qTuicCQpPaH1MwF2olIZ+AMYIaqrlHVtcAMdv0iKThffQW33QannQbHHBN1NK6QtG8Pd9wBb71lnbvXXGMdvc8+6528+a6hNf1OqroibK8EEoMMuwDLkvZbHtpqa9+FiIwRkXIRKa+urm5gePG3bRv84Ae2qtEvfhF1NK5QHXqoTeP93HM2hfe//isMHAjvvBN1ZC5TGt2Rq7ayetqODVR1vKoWq2pxUVFRup42dm64AR59FH77W1vowrkonXWWHfXffjvMmwdHHAH//u+wJvzGV7Vfpp98Ah9/DO+9B2+8YaOAyspsFs/ly6N9Dy41TRv4uFUi0llVV4TyTaIrqArolrRf19BWBZy0U/uLDXztnHfPPbb27aWX+ogdFx/NmtnZ4OefD7/8Jdx9t3X6NmliCb+usk+HDvYFUFyclXBdAzU06U8FRgHjwvWUpPYrRORRrNN2ffhiKANuTuq8PR0oyK7LZ56BK66A734X7rrL17518VNUZAcmP/4xPPCAJf3Wre2y5541b2/ZYvP2n3KK9QuUlkb9LlxtROv4+haRR7Cj9H2AVdgonKeBycD+wFLge6q6RkQEuAvrpP0SuFhVy8Pz/AC4ITztTar6QF3BFRcXa3l5ef3fVUzNnw8nnWR11BdfhDZtoo7IufSpqrJBCUuXwpNPWt+Ai4aIvK6qNf7mqjPpRymfkv5HH0HfvnZU9OqrsN9+UUfkXPqtXg1nnAGLFsEjj8CwYVFHVJh2l/T9jNws+PRTOPNMW73o+ec94bv8te++tgh7cbGd/eurdcWPJ/0M27gRhgyxI/0pU+CQQ6KOyLnMatcOpk+Hk0+Giy6CP/0p6ohcMk/6GbRtm3VuzZ5t09n26xd1RM5lR5s21qE7eLANXBg3LuqIXIIn/Qy69lqYPNkWR/n+96OOxrnsatnSVoI77zybZuSGG/xs3zho6JBNV4e77oLf/c6WPvQ58l2hatbMfuW2aWMnIn72mc0qu4cfbkbGk34GTJpkk1kNGgR/+IOPxXeFrUkTm1hwr73sQGjDBpvorUmTqCMrTJ7002jtWpunfOJEG575yCP+H9s5sAOfW2+1xP/LX9qR/n33+RF/FDzpp8mUKfCjH0F1Nfz85/Bf/2WLnDvnjIj9XWzdCr/6FbRtazN9+i/h7PKk30iffGITUz36KBx5JEybBr17Rx2Vc/E1dqzV9m+/3RL/b34TdUSFxZN+A6nayIQrrrBFUG680UbrNG8edWTOxZuILdry+edw002W+H3iwezxpN8AK1faqJwnn7QzD2fOhMMPjzoq53KHiHXufvEFXHedJf7LLos6qsLgSb8eVOHhh+GnP7X/rOPGwX/8BzT1T9G5emvSxAY9bNhgB1Ft2sCFF0YdVf7zvvMUbN1qR/Nnnw0XXAAHHwwLF9pPUk/4zjVcs2Z2AuOpp8LFF8MTT0QdUf7zpF+LbdvglVesk7ZLF5sy9uWXrfPp5Zd9Dh3n0qVlS3j6aRvmPHIkvPBC1BHlNz9OTaIKCxbYSJxJk2xZuJYt7Qh/xAhbUq5Vq6ijdC7/tGlj6/SefDKcc46twOULsWSGJ33g3Xct0T/6KHzwgZVszjjDRhYMGmQnlDjnMisxO2dpqR1ozZwJxxwTdVT5J+tJX0QGAn8AmgD3qWpG59/butVG2yxbtuPl44+3b69caaMJTj7Z5skZOhQ6dsxkVM65mhQVwd/+ZjPSDhwIt91mv7Zh+2RtydfJE7g1b24nRLZoUfN24rpNG7sU6klhWU36ItIE+BMwAFgOzBeRqar6Tjpf5+23bX3PZctsCbctW3a8v3Vr2H9/6NbNhloedRQMHw6dO6czCudcQ3TpYom/tNTm48+EPfawXxbJl7333rVtr71qv+y5Z25+cWT7SP9YoFJVlwCEBdQHA2lN+q1a2T9qSYkl9kSCT1zat8/NfyznCsWBB8J779kvctj+95p8nbytCps3w6ZN26933k5cvvjCTqhMvqxfD4sXb7+9YUPdMYrY+QV77bV9ypWd46ppG3b8lbLzduL6u9+1CRvTLdtJvwuwLOn2cuC45B1EZAwwBmD//fdv0It8+9vw0ksNjNA5Fwtt2kCvXtG89pYt9kXw+ec2ZcTuLuvX21Kou0vkydt1fSkkrnv2zMx7i11HrqqOB8aDLYwecTjOuQLUtKn16+Vj3162x+lXAd2SbncNbc4557Ig20l/PtBTRA4QkebACGBqlmNwzrmCldXyjqpuEZErgDJsyOb9qroomzE451why3pNX1WnAdOy/brOOed87h3nnCsonvSdc66AeNJ3zrkC4knfOecKiKjG9/wnEakGlqaw6z7AJxkOpzHiHJ/H1nBxjs9ja7g4x5dqbP+iqkU13RHrpJ8qESlX1eKo46hNnOPz2BouzvF5bA0X5/jSEZuXd5xzroB40nfOuQKSL0l/fNQB1CHO8XlsDRfn+Dy2hotzfI2OLS9q+s4551KTL0f6zjnnUuBJ3znnCkhOJn0ROVdEFonINhGpdfiSiAwUkfdFpFJErstifB1EZIaILA7X7WvZb6uILAyXjE4xXddnISItRGRSuH+eiHTPZDz1jO0iEalO+qx+mMXY7heR1SLydi33i4jcGWJ/U0T6xCi2k0RkfdLn9sssxtZNRGaJyDvhb/WnNewT5WeXSnyRfH4i0lJEXhORihDbr2rYp+F/r6qacxfgUOBg4EWguJZ9mgAfAgcCzYEKoFeW4rsVuC5sXwfcUst+G7IUT52fBXAZ8OewPQKYFKPYLgLuiuj/WinQB3i7lvvPAp4HBOgLzItRbCcBz0b0uXUG+oTttsAHNfy7RvnZpRJfJJ9f+DzahO1mwDyg7077NPjvNSeP9FX1XVV9v47dvlmEXVU3A4lF2LNhMDAxbE8EhmTpdWuTymeRHPPjwKkiWVk+Psp/pzqp6j+ANbvZZTDwkJq5QDsR6RyT2CKjqitU9Y2w/TnwLrZGdrIoP7tU4otE+DwSS7M3C5edR9w0+O81J5N+impahD1b/6idVHVF2F4JdKplv5YiUi4ic0VkSAbjSeWz+GYfVd0CrAeysUJoqv9Ow0IJ4HER6VbD/VGJ8v9ZKo4PZYLnReQ7UQQQSg+9sSPWZLH47HYTH0T0+YlIExFZCKwGZqhqrZ9dff9eY7cweoKI/A3Yr4a7fq6qU7Idz852F1/yDVVVEaltXOy/qGqViBwI/F1E3lLVD9Mdax54BnhEVTeJyKXYEc4pEceUC97A/o9tEJGzgKeBntkMQETaAE8AV6rqZ9l87VTUEV9kn5+qbgWOEpF2wFMicpiq1th3U1+xTfqqelojnyKji7DvLj4RWSUinVV1Rfi5urqW56gK10tE5EXsaCMTST+VzyKxz3IRaQrsDXyagVjqHZuqJsdxH9ZnEhcZ/X/WGMlJTFWnicjdIrKPqmZlMjERaYYl1IdV9ckadon0s6srvqg/v/C660RkFjAQSE76Df57zefyTpSLsE8FRoXtUcAuv0xEpL2ItAjb+wAnAu9kKJ5UPovkmIcDf9fQS5Rhdca2U513EFZ/jYupwIVhJEpfYH1SaS9SIrJfos4rIsdif+/Z+CInvO4E4F1Vva2W3SL77FKJL6rPT0SKwhE+ItIKGAC8t9NuDf97zXbPdDouwDlY/W8TsAooC+3fAqYl7XcW1iv/IVYWylZ8HYGZwGLgb0CH0F4M3Be2TwDewkarvAWMznBMu3wWwI3AoLDdEngMqAReAw7M4udVV2y/BRaFz2oWcEgWY3sEWAF8Hf7PjQZ+BPwo3C/An0Lsb1HLaLKIYrsi6XObC5yQxdhKsM7HN4GF4XJWjD67VOKL5PMDjgAWhNjeBn4Z2tPy9+rTMDjnXAHJ5/KOc865nXjSd865AuJJ3znnCognfeecKyCe9J1zroB40nfOuQLiSd855wrI/wMy6RHeYhF6/gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Find radius of furthest-out star\n", "max_rad = hy.hdf5.read_data(cantor_cat, 'Subhalo/MaxRadiusType', read_index=10)[4]\n", "print(f\"Furthest star in subhalo 10 is at r = {max_rad:.3f} Mpc\")\n", "\n", "# Load all star particles within that aperture\n", "snap_file = sim.get_snapshot_file(29) # Automatically construct the (first) snapshot file name\n", "stars = hy.ReadRegion(snap_file, 4, centre, max_rad)\n", "\n", "# Within this snapshot region, find stars that are in the subhalo\n", "ind_in_snap, in_snap = hy.crossref.find_id_indices(ids_stars, stars.ParticleIDs)\n", "\n", "# Verify that we found all of them\n", "print(f\"Located {len(in_snap)} out of {len(ids_stars)} particles in snapshot region.\")\n", "\n", "# Make a radial histogram\n", "radii = np.linalg.norm(stars.Coordinates[ind_in_snap, :] - centre[None, :], axis=1) * 1e3 # Convert radii to kpc\n", "hist, edges = np.histogram(np.log10(radii), bins=30, range=[-1, 3])\n", "midpoints = 0.5 * (edges[1:] + edges[:-1])\n", "\n", "# Make the plot\n", "import matplotlib.pyplot as plt\n", "ax = plt.gca()\n", "ax.plot(midpoints, hist, color='blue')\n", "\n", "# For reference, add stars as identified by Subfind\n", "subfind_ids_obj = hy.SplitFile(subfind_file, 'IDs', read_range=[sub.SubOffset, sub.SubOffset+sub.SubLength])\n", "ind_subfind_in_snap, subfind_in_snap = hy.crossref.find_id_indices(ids_stars, subfind_ids_obj.ParticleID)" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(107291,)" ] }, "execution_count": 69, "metadata": {}, "output_type": "execute_result" } ], "source": [ "subfind_ids.ParticleID.shape" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "25692345" ] }, "execution_count": 65, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sub.SubOffset" ] }, { "cell_type": "code", "execution_count": 80, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Prepared reading from 'IDS'...\n", "Loaded file offsets in 0.168 sec.\n", "46 \n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 80, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA25klEQVR4nO3deZyNdf/H8ddnxoyxG5oQQtlVxCQhiayVvVLd0R2pqLv1LuoupV+rVi1KUVSiUIjIVkpZRnbC2I1t7EsMYz6/P66LBjPmzMyZc52Z83k+Hucx5/qe61znfQ7zOdd8r+91fUVVMcYYExrCvA5gjDEmcKzoG2NMCLGib4wxIcSKvjHGhBAr+sYYE0LyeR3gfC644AKtWLGi1zGMMSZXWbhw4W5VjUnrsaAu+hUrViQuLs7rGMYYk6uIyKb0HrPuHWOMCSFW9I0xJoRY0TfGmBDiU9EXkUdFZIWILBeRr0UkSkQqicg8EYkXkdEiEumum99djncfr5hqO/3c9tUi0iqH3pMxxph0ZFj0RaQs8B8gVlUvA8KBrsBrwNuqWhnYB/Rwn9ID2Oe2v+2uh4jUdJ9XC2gNfCgi4f59O8YYY87H1+6dfEABEckHFAS2A82AMe7jw4EO7v327jLu481FRNz2UaqapKobgHigfrbfgTHGGJ9lWPRVNQF4A9iMU+wPAAuB/aqa7K62FSjr3i8LbHGfm+yuXzJ1exrPOU1EeolInIjEJSYmZuU9GWOMSYcv3TvROHvplYCLgEI43TM5QlWHqGqsqsbGxKR5boExAXNw60F+6/k5xw8f9zqKMX7hS/fODcAGVU1U1RPAOKARUNzt7gEoByS49xOA8gDu48WAPanb03iOyQNSTir9+8P773udxD+OH4dVdW7n76EjGfqpzTth8gZfiv5moIGIFHT75psDK4FZQBd3ne7AePf+BHcZ9/GZ6szUMgHo6o7uqQRUAeb7520YryUfS+b3qnezZcAw+vWDo0e9TpQ9qnDvvfDwnmf5vtjdvPFefk6e9DqVMdnnS5/+PJwDsn8Cy9znDAGeAh4TkXicPvuh7lOGAiXd9seAvu52VgDf4HxhTAH6qKr9GuUBx4/DXXec5ND6XbStvoH7Dr/B/Jenex0rWwY9uIYRI6DtCw1oPuQ26q3/ht8GzPQ6ljHZJsE8XWJsbKzatXeC27H9x7jz9hTGTSnI2wOTeeg/wq4CFzPvkjvosHag1/GyZHrvcTQdfCtvt5nGE5OuJ+XESbYWqsaaovW4YfdoRLxOaMz5ichCVY1N6zE7I9dk2ZHDyvJL23PvlE58NFh55Il8hEeG8959K7ht80D27/c6Yeb9+CN0+rgVIy99jkfGNEYEwiPDmf3cdFrvHclvv3md0JjssaJvsuTAAWjVWhi8rytFet3Bfff/s/vb6Z7iHD8O48Z5GDALVoyP544ux6lcuxAdFz1HRMGI0491frwi0SXDeXNgiocJjck+K/om0/au3UPvaxYxbx60+ebfNPq42xmP16sHb5X4Py7pd5tHCTNv0/JDlOzUhI/lPiZNgiJFzny8YEF4pf1c3pxYhfWT//ImpDF+YEXfZMrOnbCoXk/eWNWWCaOP0qXLueuIwGVXRrBlVyTbtwb/sfp9+6DNrUV4Puo16o18gjJl0l6vw+OXsiWsAiM/ORLYgMb4kRV947OtW6FJE+hz4h22vzGSNp0KpLtu+fefohtf8M3Y4L68UtLBJPq0Wsu6dXD75Lu4tF2tdNeNqRnDqHtn8n8/1mPnzgCGNMaPrOgbn2z+eT1fXDGQHduVT6dVoO7j1593/erVoW5dmDR8d4ASZl5KCvxc73HeX3A1X7+/h+uuy/g5jz0G+ZKO8H2/eTkf0JgcYEXfZCgxEb5rN4z79r/K7NHbadzYt+c9X2M0Py4qxabpa3M2YBb17Qv3xT/BvM4D6XRvSZ+eU7Uq/FDufm75vC1HdufyM9BMSLKibzI0fTo8emgAa0YupHabi3x+XuzDjXiFpxn3U+EcTJc1o59dzsCB0PaBirT+tkfGT0il2MtPcZNO5LNR6XdvGROs7OQsk6G+feGtt+DwYYiMzNxzmzZ1Dv6uXEnQnNR0YNN+qFiBEdVf4YFlvcmXL8OnnKNhQ9ixA9asIUvPNyYn2clZJlvaDr+V0UV6ZLrgA9xxu1L8rz9Y9eNGv+fKqo+GRTKA52jxfKMsF+wnHzrKvzc8y9xnJvo3nDE5zIq+ydCCQzVIrlg5S8/t0nwfs2nCruc/9HOqrDlyBAZ+UJC/2j5O9dtqZ3k7N3eO5M6Ib1j7xVyC+I9lY85hRd+c165d8MSRF9j6r35Zen6JyiUYcPVkeic8Q0oQnMw668GxXLXnR555OnuVOjwynBkDF3HP9pf49Vc/hTMmAKzom/Na8qdzclXtrO8UU/PhFqzaVszz69YkJUHpr97klWKv0bBR9g8w3HlvQS64AD54+YAf0hkTGFb0zXlFvP82OyhF7UsPZ3kb7drBHZFj2PrMYD8my7zhw6HhiZ85OPgrv2yvYEF494aJfDa1DOsm2aUZTO5gRd+c1/yjl/NDoa6UrJD1YZeFCkHvMt9R649POJ7kTQd4ctJJXn81hTpXRXJt13OmZs6yls814Mvw7nz8RUG/bdOYnGRF35zXl4mtGHfdu9nezqHXPuTKk3FMm+7NuM15j45i/IbLGXBfgl+Hjl5QI4ZF9w7m3e8uZscO/23XmJziy8To1URkcarbQRF5RERKiMg0EVnr/ox21xcRGSQi8SKyVETqptpWd3f9tSLSPf1XNcEg6UgyW1cezFZ//inNOhYjukQYI0dmf1uZlZICX0wqwdaiNWnZPZ2rqWXDY49BpeOrmf7Q+IxXNsZjvkyXuFpV66hqHaAe8DfwHc40iDNUtQoww10GaIMz/20VoBcwGEBESgD9gauB+kD/U18UJjhtnLiMvSeL0ebEhGxvKzISXqg/iSe+rseRXYG9SuWECfDx5jbs/vBbwvL5/4/bKlXgs4uepvnYBzi8P9nv2zfGnzL7G9AcWKeqm4D2wHC3fTjQwb3fHhihjrlAcREpA7QCpqnqXlXdB0wDWmf3DZics3TbBTzNS5Rpe6Vftte4bVH2aTFmjgxcP4imKH88MZZqlY5zWw5e3j/y3Teoo4sYNsJOzzXBLbNFvyvwtXu/lKpud+/vAEq598sCW1I9Z6vbll77GUSkl4jEiUhcYmJiJuMZf/p9S3neKfA0lZqU98v2ruhzLd3LzWTIjEv9sj1fxL07h9fWdeGj60bm6OUS6nWpxMWxpfjyy5x7DWP8weeiLyKRQDvg27MfU+cCPn4ZlqGqQ1Q1VlVjY2Ji/LFJk0WJf8RTt1YS4X66JH5YGHTtCr/8+Dd7tvztn41m4L/fN+KOklO55t3bc/y17qk5l0fj7gh495UxmZGZPf02wJ+qemr6iJ1utw3uz11uewKQetewnNuWXrsJQpqiDJp/NS8ffsiv2+3ebAsJJ0uxol/O7xLPmQO/zBbq/68l+Yvmz/HXi710H030F5ZP3JDjr2VMVmWm6N/OP107ABOAUyNwugPjU7V3c0fxNAAOuN1AU4GWIhLtHsBt6baZIJSwJYX7dTC7brzHr9ut1aocI0s+xJer6vl1u2k50rUHjxf6iHvvzfGXAqDKg60oz1ambL0sMC9oTBb4VPRFpBDQAhiXqvlVoIWIrAVucJcBJgPrgXjgE6A3gKruBV4EFri3AW6bCUJLlofzLbdSukMDv25XwoRdD7/Mp4vqsWVLxutn1ZL5SejWrdx0zR4KFcq510mteIkwLr9CPL/chDHn49OhLVU9ApQ8q20Pzmies9dVoE862xkGDMt8TBNoW6f/RRXCuOKKqn7f9u23w7DnNvDLwF38a9DVft8+wEtv5Gdq0alsGh3Yq7w9VGYMNae9Q/Kxn8kXZSN5TPCxM3JNmi4f25/p+VpTtKj/t125MvxQ6Dbqfdrb/xsH4n/bwfRv99Gnj7P3HUjVauXjWEoky38O3rmBTWizom/S9HL4c3x69Sc5tv2lvT6gzdGxzMuB+cW3/7sfq6jBIw8k+X/jGaj0aAeaM5Of/yod8Nc2xhdW9M05jhyByZtqEXbDOb13ftPqf1eRUr4it3ROYcdq/12aeNMmeHjDo0xr8ToXls/5ETtnK1cOKlaEX2fbzComOFnRN+dYM3MrbXQSdavl3HjzEiWcyyO8vOMedtdrxdEDx/2y3YEDYXnYFVw3tJtftpcVL0cP5PUJ1dAUK/wm+FjRN+c4NHoyk7iJK8vuynjlbKhTByo/3p7PjtxCj/sjsj3t4OZZ66jy0WM8eMtOyvvnJOIsKVn/Un462Zx1y496F8KYdFjRN+f4vsDttCr4K+UaVcjx12rwWkdKvvQ4X48SXn8+62fpLl4Mr3aaT7eUz3mkZ9YnfPGHcv/pRG8GMzvOrrFvgo8VfXOO+auK8Hfdxkh4YP579OsHT9y8mjsHVGXuk+MyfsJZfpl0mCZNYGLh29n5WzwXXx+4a/ukpXp1p/sqbtYhT3MYkxYr+uYMKckpNIwbRIuLVwfsNUXgxREVWFWyMQ+9V5VFi3x/7uwnJlD9pktpeuFKfv8dqjcskXNBfRQWBl8Xvpf/js75s46NySwr+uYMW+ds4vWkh2kaPjugrxtVPIrLl49iZ8xltGsH2zdlfGD3zTfhrjdrszzmer6YcZGn/fhnO3z9zbxzog87tgX25DBjMmJF35xh4d5KlGY7hbrfEvDXLl3aGdFz987X2FOrCUf3pn0gNCU5hWHtvueJJ5T6XSrQaPMoilUoHtiwGbjo/nYM4mHm/GG/Yia42P9Ic4YlSyAxrDQ1rinuyevXqQM3PlaNuCPVue+BsHNG9CQlweDGX3HPxI6812Emo0ZBVJQnUc+rbl0okf8If/0Q73UUY85gFwcxZ7hw/BD6lC5JwYKdPcvQ4NUOzCzagS+egWqXpfDMs86+yYED0LEj/DzvTsp0K0qfz5ohQbrbEhkJPxVoR4FvDsBncV7HMea0IP2VMV5ptuJ9uoafM09OwPXrB/d3TqT5cw2Z++Q4dvy5jbkVbmXF7D0MHxFGp+HtkTDxOuZ5LWvbl8ePvsQhG8RjgogVfXPagQNQ88Rifv93zl1zx1ci8PanRQgrWoS33g2nT5v11Dswi4lvrOauu7xO55sy3VowRVsxd67XSYz5h3XvmNOWLgUljJpXF/E6CuCM6Ln4r5+Ye7Vw/Dhs/XUD9RsX9jqWz665Bi6TFaz/5ii0iPU6jjGAFX2TyoEvJ/Iqv1G7xotApNdxAChdRliyxLkfHZ17Cj5A0aLwTVQ3jo4tDp/M8DqOMYDvM2cVF5ExIvKXiKwSkWtEpISITBORte7PaHddEZFBIhIvIktFpG6q7XR3118rIt3Tf0XjhZS4P+kmX3JRhQivo5whOtq55UaTbhrMv45+yokTXicxxuFrn/67wBRVrQ7UBlYBfYEZqloFmOEugzOBehX31gsYDCAiJYD+wNVAfaD/qS8KExz+L7w/dzXdEvQHSHOTCrfUZ9WxSpk6y9iYnJRh0ReRYkATYCiAqh5X1f1Ae2C4u9pwoIN7vz0wQh1zgeIiUgZoBUxT1b2qug+YBrT243sx2XDyJCxfDlfUsWP7/tSoodKJsWz8/GevoxgD+LanXwlIBD4TkUUi8qk7UXopVd3urrMDKOXeLwuknvJ6q9uWXvsZRKSXiMSJSFxiYmLm3o3Jso3T4/nyaCealljqdZQ85aKywlv5nqLc+A+8jmIM4FvRzwfUBQar6pXAEf7pygFOT4bulxkjVHWIqsaqamxMTIw/Nml8sHHeTi5jOVX9Pw96yHu/3U90Sfoq2/MFGOMPvhT9rcBWVT01m+kYnC+BnW63De7PUzNuJACpL31Vzm1Lr90EgRnHGlEr3xoqtb/C6yh5TrU2l7B9TyRr1nidxBgfir6q7gC2iEg1t6k5sBKYAJwagdMdGO/enwB0c0fxNAAOuN1AU4GWIhLtHsBt6baZILB4MdSoAfkDP61snndt/ST68TLrB9t/d+M9X8fpPwR8JSKRwHrg3zhfGN+ISA9gE3Cru+5koC0QD/ztrouq7hWRF4EF7noDVHWvX96Fyba+M1qw7souwH1eR8lzql4WySPyLvNn3IMznsEY7/hU9FV1MZDWKYXN01hXgT7pbGcYMCwT+UwA7E5I4tDxSMpeHO51lDxJwoSHbtxA3MqC3OR1GBPybHyeYclf+bmJSYT16ul1lDzr6usLsn49bNvmdRIT6qzom9OXOahd29sceVnTy/fwOd2J/2CK11FMiLOib6jz8QP8mL8DNkI251zeqCjXyWwS5m71OooJcXbBNcPiQ5dS/qKiXsfI0yIKRtDj+g3s2we3ex3GhDTb0w9xx49D391PsPDW17yOkuc1bux0pR086HUSE8qs6Ie4v1ac5MQJtf78AGhRKZ45KQ1Y/cF0r6OYEGZFP8QdHPwVicRQL2az11HyvNqtSnOCSFYsS/E6iglh1qcf4hYfvIS/wjpzd+OLvI6S5xUpU5hH682m8Da42+swJmTZnn6I+353Yz6u+zH5ouz7PxCuvRbmz03h+DHb2zfesKIfwjRF2bBov/XnB1D7C+awJSmGNV/O9zqKCVFW9EPYzsXbWbc3mtv/Hup1lJBRq31lvqcDC1YW8jqKCVFW9EPY8jWRPMWrFGt9jddRQkbMZaV4rcpQvl93uddRTIiyoh/C4jZewOs8ReV2Nb2OElIaN4bVv+4iJdn69U3gWdEPYdvmbKBq+aMUL+51ktByV/5v+GtfKdZP/svrKCYEWdEPYU9Mac7wk//yOkbIqXTHNfyX1/l9VbTXUUwIsqIforYlKP9JfovNHR/2OkrIqdC4PF+V+S9TlpTxOooJQT4VfRHZKCLLRGSxiMS5bSVEZJqIrHV/RrvtIiKDRCReRJaKSN1U2+nurr9WRLqn93om5835XRhPByp2a+J1lJAjAs0bHePA9AU2WboJuMzs6V+vqnVU9dQMWn2BGapaBZjhLgO0Aaq4t17AYHC+JID+wNVAfaD/qS8KE3hbx87jsvxrufJKr5OEpt7Jg5iUWJ9NC3Z5HcWEmOx077QHhrv3hwMdUrWPUMdcoLiIlMGZHHSaqu5V1X3ANKB1Nl7fZEOriX0YXuA+IiK8ThKaSt7XhfZ8z+w/C3sdxYQYX4u+Aj+JyEIR6eW2lVLV7e79HUAp935ZYEuq525129JrP4OI9BKROBGJS0xM9DGeyYzDh+HWY18wr9PrXkcJWZVbXsKcku2ZNa+g11FMiPH1giuNVTVBRC4EponIGWPNVFVFxC+9k6o6BBgCEBsbaz2eOWDePFiRUoNKt3idJHSFhUGXehtI+XEZ0M7rOCaE+LSnr6oJ7s9dwHc4ffI73W4b3J+nOicTgPKpnl7ObUuv3QTYthHTuZFJXGMn4nqq54kPeW/nLWyNP+Z1FBNCMiz6IlJIRIqcug+0BJYDE4BTI3C6A+Pd+xOAbu4ongbAAbcbaCrQUkSi3QO4Ld02E2DVJ73Jm1FPU6yY10lCW9RjfajNEmbPy+91FBNCfNnTLwX8JiJLgPnAJFWdArwKtBCRtcAN7jLAZGA9EA98AvQGUNW9wIvAAvc2wG0zAZScDG2PjWNk53FeRwl5NdpUZHvR6sz+VbyOYkJIhn36qroeOOfiu6q6B2ieRrsCfdLZ1jBgWOZjGn9Ztgx2HylAtbaXeh0l5IWHQ+/qM2H8Zvjobq/jmBBhZ+SGmK1DJvM4b9C4QbLXUQxwR/II+uz4H7t22pgFExhW9ENMxPQfeTR8EBdfYjNlBYOkF1+nCmv59Tfr4jGBYUU/hKhCz6Pv8b8OK7yOYlxX3HAhYQUL8MsvXicxocKKfgjZvBkSEqBe0yJeRzGuyEh4qfxHlB73oddRTIiwoh9CNnwwmSHcS5PaB7yOYlJpq5Ook/AD+/Z5ncSEAiv6IWTH3A20lGnUutqu9xJMdnwwlhuZzJw5XicxocCKfgh55UAfet2wgfDIcK+jmFSuahRJZCTWr28Cwop+iNi/3xmj36ixjRIJNgWilNElHqD0N4O8jmJCgBX9ELF28HR+1ibccMl6r6OYs4lQs+BGkrds5/Bhr8OYvM6KfohYueQEkZygdosLvY5i0rBx8I/01Vf4/Xevk5i8zop+iPh8Zxseiv2DQqXsIG4watjQuSzD7NleJzF5nRX9EHDiuDJvrtKokddJTHoKhx9lQcEmlBhl4/VNzrKiHwLWfDmf9cfK0C7mD6+jmPQUKEC+0hewfGNhjtnl9U0OsqIfAhaujGIaLah1UyWvo5jz2PTWOD472Y1587xOYvIyK/ohYPyG2jx/yReUql3a6yjmPBo1AkH5dZZdAdXkHCv6eZymKMtm77P+/Fwg+tBmEvOVJmLM115HMXmYz0VfRMJFZJGI/OAuVxKReSISLyKjRSTSbc/vLse7j1dMtY1+bvtqEWnl93djzrF5Zjxrdpege/5RXkcxGSlbljVVb2Z6fEVOnPA6jMmrMrOn/zCwKtXya8DbqloZ2Af0cNt7APvc9rfd9RCRmkBXoBbQGvhQROx6ADls7vLC/I8XKd+5vtdRTEbCw9k24FOmJ13LwoVehzF5lU9FX0TKATcCn7rLAjQDxrirDAc6uPfbu8u4jzd3128PjFLVJFXdgDOHrlWiHDZ9RRk+jP4flVte4nUU44Nrr4Vo9jJnZpLXUUwe5eue/jvAk0CKu1wS2K+qp444bQXKuvfLAlsA3McPuOufbk/jOaeJSC8RiRORuMTERN/fiUlT4vQlXNvgBGF29CZXuDD+d/ZSkgPfz/I6ismjMiwFInITsEtVA/IHp6oOUdVYVY2NiYkJxEvmWXtW7+b7jXV4NOVNr6MYX11xBeOv+j8mrKrCyZNehzF5kS/7f42AdiKyERiF063zLlBcRE5NtFoOSHDvJwDlAdzHiwF7Uren8RyTA+YuKcAtfEOhbp29jmJ8VbgwRx55hiWHL2XpUq/DmLwow6Kvqv1UtZyqVsQ5EDtTVe8EZgFd3NW6A+Pd+xPcZdzHZ6qquu1d3dE9lYAqwHy/vRNzjtkLCzEh8hYu71TF6ygmE5o0OE595vHbLBvCY/wvOz29TwGPiUg8Tp/9ULd9KFDSbX8M6AugqiuAb4CVwBSgj6raH7A56OTEybS5bAtRUV4nMZlRbv445tGATRNtV9/4nzg74cEpNjZW4+LivI6RKx3bfwyJLsYfVz1M0/mvex3HZEZiIh/c+gtvLGnB+j3FEJv3xmSSiCxU1di0HrMxHXlU3NJIrmIBJ3o+4HUUk1kxMRTs1oWN+4qxalXGqxuTGVb086g5f4SxjCuo09EuspYbNauWwK2Mtnlzjd9Z0c+jIkZ9we3lf8NGveZOF8eNYzRdWfHjZq+jmDzGin4elHJS+deSJ3igwOdeRzFZJF1v46k2S/k+rhxBfNjN5EJW9POgv1YLl+g6Eh74P6+jmKy68EIuaX85CdvDWG9z2Rs/sqKfB82ZA0coTL0b7fr5uVnL4vN5kPesX9/4lRX9PCjy0w95qMhnVK7sdRKTHRVXTOINnuD36X97HcXkIVb086DqS7/hjsITbXx3LiePPsIDt+xh/LSCJNtkWsZPrOjnMTt2QINjPzPvP195HcVkV3Q0N3UtzO7d8OuvXocxeYUV/TxmlntF3muaFfA2iPGLtgdG8ly+lxk71uskJq+wop/HFHr1WQZGPUu9el4nMf4Q9fNU7ir8HePGQUpKxusbkxEr+nmIKvy9NoF6ZRIIt4ko84ahQ5n/wQK2b4e5c70OY/ICK/p5yPLlcPvRYWz43zCvoxh/yZePm26CyEgYN87rMCYvsKKfh0yd4py62bKlx0GM/5w4QdEnevFirVGMHYudnWuyzYp+HtLo9fZ8Hd2bcuW8TmL8JiICfv+dZpdsYONGWLTI60Amt7Oin0f8/Tf8tq8WBa6wM7LynOXLqfRxP8LDsVE8Jtt8mRg9SkTmi8gSEVkhIi+47ZVEZJ6IxIvIaBGJdNvzu8vx7uMVU22rn9u+WkRa5di7CkG//AJPnnyFAk8/5nUUkwNKloSmTbEuHpNtvuzpJwHNVLU2UAdoLSINgNeAt1W1MrAP6OGu3wPY57a/7a6HiNTEmWO3FtAa+FBEbIyJn/z6/R6i8ivXXut1EuN3f/4JzZrRo+EqVq+GlSu9DmRyM18mRldVPewuRrg3BZoBY9z24UAH9357dxn38eYiIm77KFVNUtUNQDxQ3x9vwkCPzxszsUQ3Ctg5WXlPgQJw5AgtGxxExEbxmOzxqU9fRMJFZDGwC5gGrAP2q+qpK4JsBcq698sCWwDcxw/gTJx+uj2N56R+rV4iEicicYmJiZl+Q6Foy2Zl4PGHOdzmVq+jmJxQowbMm0fJtlfTqJH165vs8anoq+pJVa0DlMPZO6+eU4FUdYiqxqpqbIxN++STqT8JH3M/VR672esoJod17gxLlsC6dV4nMblVpkbvqOp+YBZwDVBcRPK5D5UDEtz7CUB5APfxYsCe1O1pPMdkQ/zXC6ha5hA1a3qdxOSYd96B2Fg6dnQWbW/fZJUvo3diRKS4e78A0AJYhVP8u7irdQfGu/cnuMu4j89UVXXbu7qjeyoBVYD5fnofISv5WDJ9Z7Xk08KP2KWU87KYGKhWjQpljhMba0XfZJ0ve/plgFkishRYAExT1R+Ap4DHRCQep89+qLv+UKCk2/4Y0BdAVVcA3wArgSlAH1U96c83E4oWxAkddRx/9/yP11FMTrrzTvjqK4iMpHNnmD8ftmzJ+GnGnE00iAf9xsbGalxcnNcxgtrzz8OAAZCY6IzlNnmcKmvjhapV4d134T/2XW/SICILVTU2rcfsjNxcLuzLEXS+bLUV/FBw/fXQowdVqsDll1sXj8kaK/q52L5NB3l63T30KW6zZIWEZs2gQQMAOnVyZtPaudPjTCbXsaKfi81YUJQKbKLgE729jmIC4dlnoVcvwBm6qQrff+9tJJP7WNHPxaZOhSPFylK3bWmvo5hASU6GpCQuuwyqVLEuHpN5VvRzKU1Rrhr1OL3r/E6+fBmvb/KAdeugcGEYMwYRZ29/1izYu9frYCY3saKfS8X/toNbDw+lddllXkcxgXLxxc5wnRo1AKfoJyfDxIke5zK5ihX9XGrSn2W4gN1UeK57xiubvCEiAl5/HerWBaBePed7wLp4TGZY0c+lpk6FS6vmo0K1KK+jmEBShU2bABBxRvH89BMcOuRxLpNrWNHPhY7tP8aTPzXn4epTvY5iAm3QIKhYEXbtApwunqQkmDTJ21gm97CinwstnLiNIikHiL3SrmIRclq2hI8+gshIABo2hNKlrYvH+M7GfeRC3y+9hHcj4tj3X6+TmICrUeP0gVyAsDDo2BGGD3fmSS5Y0MNsJlewPf1caOoUpXFjKFTI6yTGE7t2wYoVpxc7d3YK/k8/eZjJ5BpW9HOZHX9u46flZehd3sbphazu3Z2rbrqaNIESJayLx/jGundymTnTj3KM5tS7uZLXUYxXnnkGTv5zPCciAtq3d+bOPX78dHe/MWmyPf1cZsyiS3m81FdU7XSZ11GMVxo3huuuO6PpllvgwAEbxWMyZkU/F0lJTiFu6h5atnQO4JkQdfIkzJkDa9acbmrRAsqWhU8+8TCXyRV8mS6xvIjMEpGVIrJCRB5220uIyDQRWev+jHbbRUQGiUi8iCwVkbqpttXdXX+tiNippJn018g/Wb0vhnsu/MHrKMZLKSnOZZY//fR0U7580KMHTJkCmzd7mM0EPV/2F5OBx1W1JtAA6CMiNXGmQZyhqlWAGe4yQBuc+W+rAL2AweB8SQD9gauB+kD/U18UxjfTl5WiPy9Q69/1vY5ivBQR4ZyS/eijZzTfc4/zc+jQNJ5jjCvDoq+q21X1T/f+IZxJ0csC7YHh7mrDgQ7u/fbACHXMBYqLSBmgFc78untVdR8wDWjtzzeT142ZV55JVz5LTK0LvY5ivNa0KZQpc0ZThQrQurVT9JOTvYllgl+meoZFpCJwJTAPKKWq292HdgCl3PtlgdRTNm9129JrP/s1eolInIjEJSYmZiZennZw22FSfp9L6xZ2Fq4BduyAESNg//4zmnv1goQE+PFHb2KZ4Odz0ReRwsBY4BFVPZj6MXVmV/fLDOuqOkRVY1U1NiYmxh+bzBNWvjeD305ew61lfvU6igkGK1Y44/UXLDij+cYbnT8AhgzxKJcJej4VfRGJwCn4X6nqOLd5p9ttg/tzl9ueAJRP9fRyblt67cYH3+5sQvf8o6jZs6HXUUwwuOYaWLXKOaCbSkSE07c/eTJs2ZLOc01I82X0jgBDgVWq+laqhyYAp0bgdAfGp2rv5o7iaQAccLuBpgItRSTaPYDb0m0zGTh4EIaOi+ZY+9uILGxn3hici+xUrw7h4ec81KOHcwXmYcM8yGWCni97+o2Au4BmIrLYvbUFXgVaiMha4AZ3GWAysB6IBz4BegOo6l7gRWCBexvgtpkMzO45gnYHRvDfJ/zSg2byil9/hffeO6e5UiVn3P7QoWecuGsMAOJ0xwen2NhYjYuL8zqGp44ehbnFW1G0sFJvj11Ry6Ty1FNO0T94kLMnSh47Frp0cc7QbdvWo3zGMyKyUFVj03rMzusMcp9/Ds2OT+HoZ6O8jmKCTb9+zqzo+c69hFa7dlCqlB3QNeeyoh/Eko8l885rSTRoIDS6uYTXcUywKV4cotKeLjMiAv79b/jhB2cIpzGnWNEPYvMeHcWMTZfyYo+NiHidxgSlt96CkSPTfKhnT6dP/7PPApzJBDUr+kEqJQU+mlKR+cVb0ezui72OY4LVl1+meybWpZfCDTc4l+hJSQlwLhO0rOgHqYkT4cuNjTn2wVDC8tk/k0nH3LnwxRfpPnzvvbBpE0ybFsBMJqhZNQlCmqIsffQzrqh4kFtv9TqNCWoZzJjSoQPExNgBXfMPK/pBaMGQRTy74R4GNfw6rYEZxvwjIcE5G2v+/DQfjoyEu++GCRNg+/Y0VzEhxop+EHp6TF1al1zA1R/YlAMmA5GRTl/gxo3prtKzp3PVzc8/D1gqE8Ss6AeZBQtgxgxo/lQsUcXTHo5nzGkxMbBzJ+frB6xaFa6/3plVyw7oGiv6QeZw5+68GPUS99/vdRKTa/gwnrdXL9iwwdmhMKHNin4QWbnsJNu2JNOwfjJFinidxuQaEyY4u/LHj6e7SseOULKkHdA1YIcJg8hrb4QzpuBXbBoTvNdDMkEoORmSkmDPnnNm0zolf37n8vuDBjm9QaVKpbmaCQG2px8ktizYwS9fbqFXL7ggxk6/NZnQqRP8/nu6Bf+Ue+91vh+GDz/vaiaPs6IfJNb3eInlKTV4otfBjFc2JguqV4cmTeyAbqizoh8Edu6EXmueYGTTIZStUdTrOCY3evhhuOuuDFfr1Qvi4+Hnn3M+kglOVvSDwDvvwNrjFWj68R1eRzG5VYkSzvDNDHTuDNHRznfEH38EIJcJOr5MlzhMRHaJyPJUbSVEZJqIrHV/RrvtIiKDRCReRJaKSN1Uz+nurr9WROysI9eBTfup/mZPHmyzjqpVvU5jcq3+/Z0rbmYgKsrp09+zBxo2hH/9C7ZuDUA+EzR82dP/HGh9VltfYIaqVgFmuMsAbYAq7q0XMBicLwmgP3A1UB/of+qLItRNeiGOzidG0fuuQ15HMbmdqlPNM3DzzbBmDTz9NIwZA9WqwYAB8PffAchoPJdh0VfV2cDZc9m2B06NARgOdEjVPkIdc4HiIlIGaAVMU9W9qroPmMa5XyQh58gReOSHG+jeYjvVu9bxOo7J7bp1g2bNnOKfgcKF4aWXYNUqZzrF/v2dA72jRvn0dJOLZbVPv5Sqnrp80w7g1KjfssCWVOttddvSaz+HiPQSkTgRiUtMTMxivOB38vhJBrb8id274fHn7Uws4wedO0Pv3pkamlOpEnz7rXNgt0QJuP12uPZaWLgw52Iab2X75CxVVRHx276Bqg4BhoAzMbq/thtMVGFMq094/vcHqPnoHzRs2MDrSCYv6NAhy0+97jqn0A8bBs88A1dd5Vyd8+WXnRO5jh51un+OHHF+nrqlXr7mGmfiFhPcslr0d4pIGVXd7nbf7HLbE4DyqdYr57YlAE3Pav85i6+d6739Njz1cw+OtS1B97es4Bs/OnECvvvOqcDly2e8firh4c4JXLfeCi++6Jy9O3y47384FC0KkyZB48ZZyG0CRtSHDjwRqQj8oKqXucsDgT2q+qqI9AVKqOqTInIj8CDQFueg7SBVre8eyF0InBrN8ydQT1XPPlZwhtjYWI2Li8viWwtOs174hXbP16VV5yJ88w2E2aBZ40+bNzt9NgMGOLvs2bB2rVP0w8OhYEHnVqhQ2veTk+GOO2DLFhg/Hlq08NP7MVkiIgtVNTbNB1X1vDfga2A7cAKnL74HUBJn1M5aYDpO0QcQ4ANgHbAMiE21nXuAePf274xeV1WpV6+e5iV/TEzUQxTS70r10r//9jqNybPmz1dNTg74y+7YoVq7tmpkpOq4cQF/eZMKEKfp1FWf9vS9kpf29FevdsZF31xgOm/NupISVUp6HckYv9u3zxkNtGCBM2nLv/7ldaLQdL49fetcCIDE5Tt5+vo/CA+HZ3+5wQq+yXmjRzsHdgO8Uxcd7UzCft11zlUhBg8O6MsbH1jRz2FHjsDC6x5lyPabmDz6kI1uMIFx+DDs3g17z3vYLEcULuwc0L35ZmcE6euvBzyCOQ8r+jkoORm6doVu+wbx1/+NJfZ6G49vAuSee+C335yZUzwQFQVjxzr//596Cv73PzvpK1jYJCo5RFOUoTd9z6Sp7Xnv/Qto1Kep15FMKDk1heKRI84wzuLFAx4hIgK+/PKfs38PHnQuLmgj1rxlH38O+abXdO6b2omvbvyaPn28TmNC0pEjzvDNV17xLEJ4uDNF42OPwXvvQY8ecPKkZ3EMtqefI4YNgx5Db2Bz4/E8/v1NXscxoapQIejXDxp4ewKgCLzxhnPy1vPPO4cbRo50/hIwgWdF3492rUjkrzaP8vKW57n++so8NLkdYfYJGy89+qjXCQCn8PfvD0WKwOOPO38BfPWV89MElnXv+IGq03fZrPFxqm6Zzhu3L+Snn5yDWcZ4bvt2ePfdoDiS+thjzmie0aPhvvts2kYvWNHPpu1xCQy97C3uuguKVC/L/oXr6TDyNvLZHr4JFlOmOHv8y5Z5nQSA//4Xnn0Whg51vgSC4LsopFhpyiJV54zDTQ98wZNJA+Dpjvx7QCXCwwt6Hc2YM526XnLlyl4nOe2FF5zRPO++6/T1DxjgdaLQYXv6WZDwx2Z6N1rCPffA7KseZ+f05fR8qZL1T5rgFBUVVAUfnD7+t992RvO8+CIMHOh1otBhe/qZkJICQz5WGj94Ez3Jz+Xvz+f+ByIIC7vE62jGnJ8qPPCAc3H8F17wOg3gFP6PP3ZG8zz5pHOQ9/77vU6V91nR90HSwSQWvzqF52c3Y8qcIvS56lP6vnUhvRuL19GM8Y0IHDvmzIYSRMLD4YsvnFMKevd2TuSyi7TlLCv66Ug+lszPU47x1fjCbP92AVOOdKB8kZF88snt9OhR//QJj8bkGp99RjD+x42IcKZsvPFGZ7auQoWgY0evU+Vd1qefSkoKzJ4Nj9x7hL0FyzG741uMGwcXdWlI3Ms/8cGOLvTsGZS/N8Zk7NR/3PXrg27ITFSUM/nKVVfBbbfB1KleJ8q7Qn5PXxXi4mDPwwPYuOwQDxweSIEChWhU8wFa3tqYp5+EqKgwwKYCMnnA1KnQujXMnAnXX+91mjMULgyTJ0OzZs6e/tSpzqAj418BL/oi0hp4FwgHPlXVV3Py9f7+25nCbe+vK0iKW8bsi7qyZQvUn/0Gl2+dTIuwmRw+DB+G7aT6RQf5+hO46SYoXLh/TsYyxhvXXedc/axWLWf55ZedtkaNvM3lio52in2TJk53zyuv/HOSo+o/f6Ccup/6D5bISMif/59besuFCzvXn4uKCs2/2gNa9EUkHGc6xRY4Uy8uEJEJqrrSn6+zZAlMvvEDem57gTK6jZPkYwCjeJqXuYEuXFAqH5WiipISfQE9Oiu16wgd2r9PdIkQ/B9gQktUFDz9tHP/yBFnvCQ4RT8pCTp1gocfhpYtPYt44YUwfbrzXfTggzn3OpGRUKyY8wWQ+pa6rWjR898KFsyFXxzpzaOYEzfgGmBqquV+QL/01s/qHLkbNqg+W3+K/lLzfn39uUP6xReqc8Zu102z1umxv09maZvG5ElJSaqHDjn3N25UveKKfya4Xb1aNSZGdfp0Z3nBAtVKlVR//dVZ/uMP1apVnTl5VZ32mjVVFy92lmfMcJZXrnSWp0xxlteudZYnTnSWN250lseOdZa3bVNV1eNfjNKkKjV1y6JE3bxZdc87I/R41ZqasHK/JiSo7n1zqB6vWlO3rzui27ap7n7pIz1WuaYu+/O4xsWprntkkB6uWEunTFEdP151Ubc3dU+FOvrRR6qvvqo6uemruv7C+tq1q2qbNqofl39R4wo00jJlVAsUUH2e53QmTU//TfEyfXUKLU8vD+Rx/YEbtVgx1QsvVP2kwEP6Q/6OWrq0aunSqp8XvF/HR92qF12kWras6siCPXRswTu1XDln+dsCd+moAndrmTLO+uOiuuqIAvdqqVLO9nr1yvo/K+eZIzfQ3TtlgS2plrcCV6deQUR6Ab0ALr744iy9SMWKMGBeK6AVTU63ls7StozJ0yIjnRtAhQrOn8mnqEL79hAT4ywXLer8RXDq2vyFC8OVVzo/wRl2U7MmFCjgLBcp4iyf6p8pWtRZzp/fWS5WzFk+9frFizvL7jVMIi6Mhto1KVchHKKByiXgippcVD4cCgOXloQralL6ojCIAqpdAHVqctnl4vRhxF8ICTVp1cp9P4dKwYka3Hefu1yqFMyqztfD3eVPS8Mf1dg21FlMfr8MyQuqsv555+zhwp9dRHj8QUbe6SxX+aEsBXaeoFsDZ8qC6EXliTgaSbsGzkdXZNHFRBwvRpurnOWoxRXIl5xES3fm2vx/VkTDwrmxrvPXQtSCSki+QnSo5zwem+YMt9kX0InRRaQL0FpVe7rLdwFXq2qaf8TlpYnRjTEmUIJpYvQEoHyq5XJumzHGmAAIdNFfAFQRkUoiEgl0BSYEOIMxxoSsgPbpq2qyiDwITMUZsjlMVVcEMoMxxoSygI/TV9XJwORAv64xxhi7DIMxxoQUK/rGGBNCrOgbY0wIsaJvjDEhJKAnZ2WWiCQCm3xY9QJgdw7HyY5gzmfZsi6Y81m2rAvmfL5mq6CqMWk9ENRF31ciEpfe2WfBIJjzWbasC+Z8li3rgjmfP7JZ944xxoQQK/rGGBNC8krRH+J1gAwEcz7LlnXBnM+yZV0w58t2tjzRp2+MMcY3eWVP3xhjjA+s6BtjTAjJlUVfRG4RkRUikiIi6Q5fEpHWIrJaROJFpG8A85UQkWkistb9GZ3OeidFZLF7y9FLTGf0WYhIfhEZ7T4+T0Qq5mSeTGa7W0QSU31WPQOYbZiI7BKR5ek8LiIyyM2+VETqBlG2piJyINXn9lwAs5UXkVkistL9XX04jXW8/Ox8yefJ5yciUSIyX0SWuNleSGOdrP++pjePYjDfgBpANeBnIDaddcKBdcAlQCSwBKgZoHyvA33d+32B19JZ73CA8mT4WQC9gY/c+12B0UGU7W7gfY/+rzUB6gLL03m8LfAjIEADYF4QZWsK/ODR51YGqOveLwKsSePf1cvPzpd8nnx+7udR2L0fAcwDGpy1TpZ/X3Plnr6qrlLV1RmsVh+IV9X1qnocGAW0z/l04L7OqZk3hwMdAvS66fHls0ideQzQXEQkSLJ5RlVnA3vPs0p7YIQ65gLFRaRMkGTzjKpuV9U/3fuHgFU4c2Sn5uVn50s+T7ifx2F3McK9nT3iJsu/r7my6PsorUnYA/WPWkpVt7v3dwCl0lkvSkTiRGSuiHTIwTy+fBan11HVZOAAUDIHM2UmG0BntwtgjIiUT+Nxr3j5/8wX17jdBD+KSC0vArhdD1fi7LGmFhSf3XnygUefn4iEi8hiYBcwTVXT/ewy+/sa8ElUfCUi04HSaTz0jKqOD3Ses50vX+oFVVURSW9cbAVVTRCRS4CZIrJMVdf5O2seMBH4WlWTROQ+nD2cZh5nyg3+xPk/dlhE2gLfA1UCGUBECgNjgUdU9WAgX9sXGeTz7PNT1ZNAHREpDnwnIpepaprHbjIraIu+qt6QzU3k6CTs58snIjtFpIyqbnf/XN2VzjYS3J/rReRnnL2NnCj6vnwWp9bZKiL5gGLAnhzIkulsqpo6x6c4x0yCRY7+P8uO1EVMVSeLyIcicoGqBuRiYiISgVNQv1LVcWms4ulnl1E+rz8/93X3i8gsoDWQuuhn+fc1L3fveDkJ+wSgu3u/O3DOXyYiEi0i+d37FwCNgJU5lMeXzyJ15i7ATHWPEuWwDLOd1c/bDqf/NVhMALq5I1EaAAdSde15SkRKn+rnFZH6OL/vgfgix33docAqVX0rndU8++x8yefV5yciMe4ePiJSAGgB/HXWaln/fQ30kWl/3ICOOP1/ScBOYKrbfhEwOdV6bXGOyq/D6RYKVL6SwAxgLTAdKOG2xwKfuvcbAstwRqssA3rkcKZzPgtgANDOvR8FfAvEA/OBSwL4eWWU7RVghftZzQKqBzDb18B24IT7f64HcD9wv/u4AB+42ZeRzmgyj7I9mOpzmws0DGC2xjgHH5cCi91b2yD67HzJ58nnB1wBLHKzLQeec9v98vtql2EwxpgQkpe7d4wxxpzFir4xxoQQK/rGGBNCrOgbY0wIsaJvjDEhxIq+McaEECv6xhgTQv4fAF2TM2QidF8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# For reference, add stars as identified by Subfind\n", "subfind_ids_obj = hy.SplitFile(subfind_file, 'IDs', read_range=[sub.SubOffset, sub.SubOffset+sub.SubLength])\n", "ind_subfind_in_snap, subfind_in_snap = hy.crossref.find_id_indices(subfind_ids_obj.ParticleID, stars.ParticleIDs)\n", "\n", "# With Subfind, we cannot pick out the star IDs separately, so some of these will not have been\n", "# matched (all those that are not stars). We therefore need to extract specifically those that have\n", "# been matched, with the second `subfind_in_snap` index.\n", "radii_subfind = np.linalg.norm(stars.Coordinates[ind_subfind_in_snap[subfind_in_snap], :]\n", " - centre[None, :], axis=1) * 1e3 # Convert radii to kpc\n", "hist_subfind, edges = np.histogram(np.log10(radii_subfind), bins=30, range=[-1, 3])\n", "ax = plt.gca()\n", "ax.plot(midpoints, hist, color='blue')\n", "ax.plot(midpoints, hist_subfind, color='red', linestyle=':')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That looks good! In the centre, Subfind and Cantor give near-identical profiles. The difference is that Cantor (blue) has a more extended tail outside ~20 kpc, as expected." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Example 3: Select only star particles within 30 pkpc\n", "\n", "The start point is the same as in example 2 (the innermost star particle), but we now retrieve the end point from the `OffsetTypeAperture` list in `Extra`:" ] }, { "cell_type": "code", "execution_count": 83, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Retrieved 79780 IDs.\n", "First 10 IDs are: [348327033 151511633 50898603 24679775 209739697 348974305 271038371\n", " 275338903 44062981 255913957]\n" ] } ], "source": [ "offset = hy.hdf5.read_data(cantor_cat, \"Subhalo/OffsetType\", read_index=10)[4] # As before\n", "end = hy.hdf5.read_data(cantor_cat, \"Subhalo/Extra/OffsetTypeApertures\", read_index=extraID)[4, 2]\n", "\n", "# Note: we again select the \"end\" index directly from the `OffsetTypeApertures` list.\n", "# 30 pkpc corresponds to index 2 in the third dimension (see above). As always, we truncate\n", "# the first dimension by loading only the entry for \"our\" subhalo with `read_index`.\n", "\n", "ids_stars_30kpc = hy.hdf5.read_data(cantor_ids, 'IDs', read_range=[offset, end])\n", "print(f\"Retrieved {ids_stars_30kpc.shape[0]} IDs.\")\n", "print(\"First 10 IDs are:\", ids_stars_30kpc[:10])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, we got (slightly) fewer IDs this time, and the first 10 IDs are the same as before (since these are the ten most central star particles). Let's again compare histograms:" ] }, { "cell_type": "code", "execution_count": 84, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Located 79780 out of 79780 particles within 30 kpc in snapshot region.\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 84, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3gElEQVR4nO3dd3hUZdrH8e+dRgIhCV0gkVAiRaWGSEeko1JEBZWmICrqrrKubXdl7b2+igsCChZEigpKFUEEpYTeJYCU0AIkoSck87x/zAEDJGQSZuZMMvfnuubKmee0XwZy5+Sc55xHjDEopZTyDwF2B1BKKeU9WvSVUsqPaNFXSik/okVfKaX8iBZ9pZTyI0F2B7ic8uXLm9jYWLtjKKVUkbJy5crDxpgKuc3z6aIfGxtLYmKi3TGUUqpIEZFdec3T0ztKKeVHtOgrpZQf0aKvlFJ+xKWiLyKPi8hGEdkgIhNFJFREqovIMhFJEpFJIhJiLVvCep9kzY/NsZ1nrPatItLZQ9+TUkqpPORb9EWkKvA3IN4Ycx0QCPQFXgfeNcbUAlKBwdYqg4FUq/1dazlEpJ613rVAF2CkiAS699tRSil1Oa6e3gkCwkQkCCgJ7AduAqZY88cDPa3pHtZ7rPntRUSs9q+NMRnGmJ1AEpBwxd+BUkopl+Vb9I0xycBbwG6cxT4dWAmkGWOyrMX2AlWt6arAHmvdLGv5cjnbc1nnPBEZKiKJIpKYkpJSmO9JKaVUHlw5vVMG51F6daAKUArn6RmPMMaMNsbEG2PiK1TI9d4Cpbxmb8oxhnz4GSdOZ9odRSm3cOX0TgdgpzEmxRhzFpgGtASirNM9ANFAsjWdDMQAWPMjgSM523NZRxUDWdkO2v53BHe8+aHdUdwiMxPqv9KLsUfu5b6PPrE7jlJu4UrR3w00E5GS1rn59sAmYAFwu7XMQOB7a3q69R5r/s/GOVLLdKCv1bunOhAHLHfPt6HsdiYzi9pPDWKRvMCUtGc4euy03ZGuiDFw//2QOuVlAL498A6ZZ7NtTqXUlXPlnP4ynBdkVwHrrXVGA08Bw0UkCec5+7HWKmOBclb7cOBpazsbgW9w/sKYDTxsjNGfomIgMxPu7pfNjgMpVE3vDSEneGnyD3bHuiJDR6xiwgR4/v5mPBEzhazSO3hmwrd2x1LqiokvD5cYHx9v9Nk7vi3txBnuvOcU86aX5c23s/jbo0LYv2K4+syt7PxglN3xCmXg+2OYkHY/HVK+Y+7/9eBsVjalnq5NqTO1Sf3wR0TsTqjU5YnISmNMfG7z9I5cVWgHU09Q4z83M69SFz76OIsnhgcREhzIfdnLSB71P9LS7E5YcC9MnMWEow9SLq0z37/VDREICQ7kmWozSB89lcWL7U6o1JXRoq8KZfehdOJe6ExqxEIeaPgowx7864GtD9wVw9lMYdo0GwMWwpc/r2LEhjsIO1afdf+eTMnQ4PPznh5Sl3KRobz5lsPGhEpdOS36qsC27T1C3Vfbc7z0coZf/TX/G9b/gvlNmkDZni/xzKo+NiUsuA3bjjNg9i0EZpbjt0d+pEq50hfML1kSejy8lBmxccxcvsWmlEpdOS36qkAOHoQmLwzhVPgGnqv9LW8PvuOSZUTgugZnOVR+Mmu277chZcGkpsKdPUoT+uvrTOs9i4Y1K+e63PB7a0DpffxjyjteTqiU+2jRVy7buxfatIGzM97jrYazeP6eW/Jc9qlud4EYnp8yyYsJC+7YyQw69FvN9u0w87X+dG9WL89lr42tSN3MgWwpMYENOw96MaVS7qNFX7lk4dodXPvY0+w/4GDe5Gr8o3e7yy7fLaEOYWmNmHfgKy8lLLisbAfX/Wcgqxq15N0xybRtm/867945HAIzeeiz4nEDmvI/WvRVvjbvTqH9F204XusTJsz4k1atXFvvpop3czJqBfNXJ3k2YCG1eO5p9kROomvYCIb1v+QxULnqHH8NldN7sCRjJIdST3o4oVLup0Vf5WvkrJ9whCfzbvNp9GxTw+X1RtzWFxIf5McZvvcE7bvf+R8rQt7k2lMP8cMzTxZo3Vc6j8B8N46vvwjzUDqlPEdvzlL5avbvp1km73L82eOEh4UUaN0bb3Re/N20CZ+5qWnXwTRi369G2VM3kPzGTEJDgvJf6SItWsCBA/DHHxBU8NWV8ii9OUtdke0n1xB6ol6BCz7AXXcZtpz4nVm//+n+YIU07pMQWPgcI3u9WaiCD/Do8NPsjP0P/5oww83plPIsLfoqX6eTbqCe485Crdv+llS4tw3//XGkm1MVzsmT8NF7JelW5h/0adug0Nvp3SOE4EZf89H6F3E4fPevZaUupkVfXdahQ3Dyh+fpF/tModavVbUsFY93YVXmRLKy7b+b9cGPvuJI9HieffbKCnVIcCC9qwznZNQKPvrhVzelU8rztOiry1q26hQEZNGg8AfF3FnnbrLD9/Lxj/Y+uObYyQy+OvxPom78jJYtr/wCw4dDBiGny/PygrfckE4p79Ciry5r1OqR8Gw41eukF3ob/7mzO2SW5OPF9vbZH/bJeByl9vFUy2fdsr1ykWG0CX2Eg1Ez+GHZZrdsUylP06KvLmvT4bUEnClP9SqRhd5GxTKlqHamB1uz5pGRYc/57zOZWUza+zql0pryZO8Obtvux/cNI2DznXz2qY90TVIqH1r01WXtN2spl3UF53YsL7Z6F8cHm/jpJ3uK4/Bxk8gqvYNHGz5LQID7MtS9ugJDy0xixqd1OHDAbZtVymNcGRi9toisyfE6JiKPiUhZEZknItusr2Ws5UVEPhCRJBFZJyKNc2xroLX8NhEZmPdelS84djKDM+GbqRV+5UW/z82VKBtZgq9sOMPjcMCPU8oSsed2XuzX3e3bHz4cMiO28uiH3+e/sFI2c2W4xK3GmIbGmIZAE+AU8C3OYRDnG2PigPnWe4CuOMe/jQOGAh8DiEhZYARwA5AAjDj3i0L5ppkrNkFgFglXN7zibYWEQEK/H5kY0cTrjy+YPh12z+/KyHaTCQp0/x+3cXFQpd+zTM2+lwNHT7h9+0q5U0F/AtoD240xu4AewHirfTzQ05ruAUwwTkuBKBGpDHQG5hljjhpjUoF5QJcr/QaU5+zbXh7mv0yPJs3dsr1u7UtjrlrFi99Md8v2XOFwGB7/fCyxtY/Rx4OP9x/R8QlMaCoPjxnnuZ0o5QYFLfp9gYnWdCVjzLmHpR8AKlnTVYE9OdbZa7Xl1X4BERkqIokikpiSklLAeMqd9myIISzxWdo0iHHL9h66uRWBJ6KZvGVi/gu7yRtTf+LP+kNo/eAkjz4uYWi35pRMi2fO/i88txOl3MDloi8iIUB3YPLF84zzAT5u6ZZhjBltjIk3xsRXqFDBHZtUhfTr9kTqND5MoJuelxYUGECjkL4cLD2LbXuPuGej+XhtySsEnKjKh/cP8Pi+GoR35GTp1fr0TeXTCnKk3xVYZYw5N3rEQeu0DdbXQ1Z7MpDz0DDaasurXfkgh8Owql5nTrdwT5/2c4Z3vBsCs3h+ylS3bjc3//vxN9LLLKR7+SeIKFXC4/vrUq81ZIcw7detHt+XUoVVkKJ/F3+d2gGYDpzrgTMQ+D5H+wCrF08zIN06DTQH6CQiZawLuJ2sNuWDVm5LxoQd5fqKV95zJ6c+bRoSte0h1i+Mc+t2c/PcT68gp8szauj9Ht8XwEOdOsDraRxc3Tj/hZWyiUtFX0RKAR2BaTmaXwM6isg2oIP1HmAmsANIAj4BhgEYY44CLwIrrNcLVpvyQT8krgHgxjruLfoBAcLw2iNZP70de/bkv3xhLV+ZQcrhs7Qv9RgVy5Ty3I5yqFAumPrXBrPY3qdNKHVZLhV9Y8xJY0w5Y0x6jrYjxpj2xpg4Y0yHcwXc6rXzsDGmpjHmemNMYo51xhljalmvT93/7Sh3+W3HWgC631Df7du+6y4wkTt586tlbt/2OW+9XoKI6XOY9EjhHhRXWJXbT+Hn2Facyczy6n6VcpXekatytSV1LUHHqxNdIcLt265VC0oN6sOYA8Pcvm2AX1bvZfLc3Tz8MJQt493/4tfUycIRvYTJv6716n6VcpUWfZWroMUjuOHgJx7bfser7uZ01CrGznH/0f6QL/8Dw67n/oe934umf1vnAMLTEvVxy8o3adFXlzh5EnYlXkuHmu09to+3+g8g8EQMD/x0O2u2789/BRct2biLpJJf0MAMonpV75zLz6lp7WiCjsey/ICe2Fe+SYu+usTc5X9i6o+n1nVpHttHzSpl+fKW6WSHHKXVh71IPX7GLdt96PM3AWHUgCfcsr3CiDGtOBCyWEfUUj5Ji766xNS1c6HXIGKu8Wznqj5tG/JU3Jec3Niahx4IxlxhjVy4dgfrg8dwzekB3FDXPXcRF8ZNV3fDsbM167boTVrK92jRV5dYu38tZETQ6rpYj+/rtUE9efnGN5k0MZD/vnyq0NtZswZ6PbIMySrJqH7uvaGsoIZ3vAsmT2bV0nBbcyiVGy366hK7MtYScao+gQHe+e/xzDNw66CtvJB6DU9+Oi3/FS7y3ZyjtGkD4X/exeI7k7ixQQ0PpHRdnTpQtiwsWHLc1hxK5UaLvrpAVraD4yXXUS3UvTdlXY4ITPi/aoQ7YngzqT8TF652ed1HR02k18LqlG/4O7/9Bi0alfVgUtcEBED43UOYFNXE7ihKXUKLvrrAkg27oMRxGl7lvaIPEBUeypJHvyUwsxz9f+zuUo+eW199mw8P3E3k6YYsmFyXGPtO41/iusq1ORuxjfU7Dua/sFJepEVfXeDojurw1n4GNr3D6/uuX+OqC3r0HD12OtflsrIdNHl2OD9kPkHV9Nv586U5VKsU5d2w+ejVuDUAExYusTmJUhfSoq8usHYtBJy6iuaNomzZf5+2DXky7gtOHirHA8POXtKjJyMDWj30JatKvEv904+y4/WviQoPtSXr5fRt2xjOhjF/m96kpXyLFn11gSn73qZC57GULGlfhtcH9eKlej8w5csIXnrZcb49PR26doVlY+5hQMh3rH7lfUKC3fSwfzcLDwsh8sQNbD2jN2kp36JFX11ga8SHhNSZZ3cMnn1W6D0ghed2teDJT6exats+Yp7uyqJ1u5gwPoDxz/QgIEDsjnlZnSMe5/TPj3NcO/EoH6JFX523+1A6WaX/pG4Z717EzY0IjPmoNKXC4c2k/iSMbsbxsot5/X9/0r+/3elcc1/L7ph1d7N0qd1JlPqLFn113vRl6wBoUdP+og/OHj2/PfodgZnlMAGZfNH+F/5xe1u7Y7mseXOQShv5ZnFi/gsr5SUeHCpaFTW/bHE+DvjmeN8o+uDs0bNtuDNX9cplbE5TMBEREHrXAKamR/EJ8+2OoxTg+shZUSIyRUS2iMhmEWkuImVFZJ6IbLO+lrGWFRH5QESSRGSdiDTOsZ2B1vLbRGRg3ntUdth+4DByogqNa1WxO8oFqlcuU+QK/jnXlGhFaqmlnDpz1u4oSgGun955H5htjKkDNAA2A08D840xccB86z04B1CPs15DgY8BRKQsMAK4AUgARpz7RaF8Q9Cv/+XGVXt8/gJpUXJTXCsIOcU3i9bYHUUpwIWiLyKRQBtgLIAxJtMYkwb0AMZbi40HelrTPYAJ1rCJS4EoEakMdAbmGWOOGmNSgXlAFzd+L+oKZGfDhg3QsIFe5nGnc4OqfLtK++sr3+DKT3h1IAX4VERWi8gYa6D0SsaYc/fKHwAqWdNVgZxDXu+12vJqv4CIDBWRRBFJTElJKdh3owpt7qo/ON2nLaXrLLc7SrHSqFZlgo7VZMVBLfrKN7hS9IOAxsDHxphGwEn+OpUDOAdDB9wyYoQxZrQxJt4YE1+hQgV3bFK5YPaaVRC7iHq1Q+yOUux0PT6VjMnjrni8AKXcwZWivxfYa4w5N5jpFJy/BA5ap22wvh6y5icDOR99FW215dWufMCK3WshO4iu8XXtjlLsdE9owNHkMvzxh91JlHKh6BtjDgB7RKS21dQe2ARMB871wBkIfG9NTwcGWL14mgHp1mmgOUAnESljXcDtZLUpH5B0Yi2hJ+oSUaqE3VGKnYTmGdD6FT6ep//dlf1c7af/KPCliIQAO4B7cf7C+EZEBgO7gDutZWcC3YAk4JS1LMaYoyLyIrDCWu4FY4xnx+NTLjsSvJZYh+cGQvdn19UNQZq/xw9/duM9OtsdR/k5l4q+MWYNEJ/LrEuqhHV+/+E8tjMOGFeAfMoLkg9m4EhuSJPr29gdpVgKCBCuymzFrgB9+Jqyn/bPU2zZUAK++pGh8UPsjlJsNa3UiqyI7axOyn9wGKU8SYu+YvUa5+OLG/jO0xeKnZ6Nnf31P/9Fj/aVvbToK/63fzDBD7RGe8h6Tp82jeBUORK3aoc1ZS994JoiOXs1kWGV7Y5RrJUMDaZd4kHSUn1z0BflP/RI38+dOJ3JmdKbqFlKz+14WutWgaxdC8eO2Z1E+TMt+n5uVuJmCDxL05iGdkcp9qo3ScJxXzNGztbHLCv7aNH3c/PWOZ9V30mv4npc55ZXQZVEftz4i91RlB/Tou/nju2qQcDqoXRsHGd3lGKvcrlwSh5rxPp0ffiaso8WfT93eFUrGiePIjREr+l7Q+3QVqSHL+PE6Uy7oyg/pUXfjzkchlU7d1C/gcPuKH6jfVwrCD7NpF9W2x1F+Skt+n5szfb9pA6oybHaH9sdxW8MbNcKtvRg7Rr90VP20L/p/dgPic6LuG1qX29zEv9xXfVKxK38jl2n7E6i/JUebvix37Y7i373G+rbnMS/tGoFv646RFa2nlZT3qdF349tTl1L4PFqVKsUZXcUv1Ki8TekDqnEzBVb7I6i/JAWfT+2T5ZRydHE7hh+5/aWjQH4+rdFNidR/kiLvp9KTjZk/fgOvSr/3e4ofqddg5oEnKzMb8la9JX3uVT0ReRPEVkvImtEJNFqKysi80Rkm/W1jNUuIvKBiCSJyDoRaZxjOwOt5beJyMC89qc877ffBLb0ZEBbHTjF2wIChOistuwJWITDoaOlK+8qyJF+O2NMQ2PMuRG0ngbmG2PigPnWe4CuQJz1Ggp8DM5fEsAI4AYgARhx7heF8r6Jy+ZTosZyGjWyO4l/alG1DY7wZBau22F3FOVnruT0Tg9gvDU9HuiZo32CcVoKRIlIZaAzMM8Yc9QYkwrMA7pcwf7VFZiT/SRhtz5NcLDdSfzTfa27wfdj2LSyrN1RlJ9xtegbYK6IrBSRoVZbJWPMubHfDgCVrOmqwJ4c6+612vJqv4CIDBWRRBFJTElJcTGeKogDR09wKmIt10W2tDuK32ofX41yuwezcrH+sau8y9Wi38oY0xjnqZuHReSCE8HWYOhuOTlpjBltjIk3xsRX0KGcPOLzBcsgIJsudVvZHcVvBQRAfPs9zNzzld1RlJ9xqegbY5Ktr4eAb3Gekz9onbbB+nrIWjwZiMmxerTVlle78rKZGxaDEfrf1MzuKH4trPFUDrW+hxVb99odRfmRfIu+iJQSkdLnpoFOwAZgOnCuB85A4HtrejowwOrF0wxIt04DzQE6iUgZ6wJuJ6tNedn61KWEpl/P1RUj7Y7i125v6vyD+bMF2nVTeY8rR/qVgMUishZYDvxojJkNvAZ0FJFtQAfrPcBMYAeQBHwCDAMwxhwFXgRWWK8XrDblRVlZcGb8NHpnT7M7it+7o3UDyIhgwQ4t+sp78n3gmjFmB3DJsErGmCNA+1zaDfBwHtsaB4wreEzlLuvXw8m0MLo1q2l3FL8XEhxIhdOt2C46kpbyHr0j18+M/Plb6PA0N7Q4a3cUBTQp34bMyC1s2X3E7ijKT2jR9zNzkicR2PBLasZqB31f8FjrIfDmQTauKGd3FOUntOj7EYfDkBy4mKoO7arpK9o1K0dJU5Ff9AyP8hIt+n7k9827cYQnk3CV3pTlK0JCoMatU5iYqg++U96hRd+PfPXrYgBuT9AjfV9Sod4mDtf8P3buT7U7ivIDWvT9yPptaUh6NXq10OERfUmPBm1BDJ/OX2J3FOUHtOj7kfSfHqbDxp2EBAfaHUXl0L9dAmSFMHuz9tdXnqdF30+kpTn76LdqKXZHURcpGxFGxIkENp3Sq7nK87To+4mPZs/BPHQt1RpvtTuKykWD0h04mV6C9GPZdkdRxZwWfT8xa9OvUH4rXVtG2x1F5eLfrUbAp4tYtlRPvSnP0qLvJzYcX0zJY42oWKaU3VFULlq0gMBA+GWRDp+oPEuLvh84eSaT9FLLuSZUu2r6qvBwKDdgGCOPdbI7iirmtOj7gUmLVkPwaW6qpTdl+bJqlcNJi1xE6vHTdkdRxZgWfT+waW0orO3HgBv1SN+XdanbBoIymfDzcrujqGJMi74f2Lm0ATXWfk6DmlfZHUVdxr3tW4IRvl+rXTeV52jRL+YcDsMv63bQoqVeIPR11SuXITS9PmtT9SYt5TkuF30RCRSR1SLyg/W+uogsE5EkEZkkIiFWewnrfZI1PzbHNp6x2reKSGe3fzfqEj+vSeJIv5oENP7M7ijKBc2DH+L4yls4q8MdKA8pyJH+34HNOd6/DrxrjKkFpAKDrfbBQKrV/q61HCJSD+gLXAt0AUaKiHZK9rCvf3M+z+W2hASbkyhXPHzDA5z99TFWrrQ7iSquXCr6IhIN3AyMsd4LcBMwxVpkPNDTmu5hvcea395avgfwtTEmwxizE+cYulqJPOzX3YuRM2W4+Ya6dkdRLmjdGgg7yrcLt9sdRRVTrh7pvwc8CTis9+WANGNMlvV+L1DVmq4K7AGw5qdby59vz2Wd80RkqIgkikhiSkqK69+JytWfWUuocKYFQYF6+aYoqFgRQh5szacH/2Z3FFVM5VsJROQW4JAxxit/cBpjRhtj4o0x8RUqVPDGLoutrXsOkxm5hUbltKtmUVIzqDUpYYvJPKvP4VHu58rhX0ugu4j8CXyN87TO+0CUiARZy0QDydZ0MhADYM2PBI7kbM9lHeUBaxPD4JtvGNC0t91RVAG0q9EGShxj6pJ1dkdRxVC+Rd8Y84wxJtoYE4vzQuzPxph7gAXA7dZiA4Hvrenp1nus+T8bY4zV3tfq3VMdiAP0LhQPWrm0FCFJd3Bb2zi7o6gCGNSuDQDfLNP++sr9ruRE71PAcBFJwnnOfqzVPhYoZ7UPB54GMMZsBL4BNgGzgYeNMfr3qwd9mzSRum02EhpqdxJVEE1rRxN0vAbLDmh/feV+Qfkv8hdjzEJgoTW9g1x63xhjzgB35LH+y8DLBQ2pCi7txBm21RtEU8ffgTfsjqMKqP3JT1j2UxXMOyA67o1yI+3SUUx9tTARgjLpUFsfslYU9Um4ibSkOmzenP+yShWEFv1iasZa501Z/W9sYXMSVRitWmdDgwl8Mm+h3VFUMaNFv5hadXgxIcdqU/dq7fZaFNWsEUBAx2eZtnuU3VFUMaNFvxjKzjakhCwnNkD75xdVAQFCjKMNewMX4XDow/KU+2jRL4a2bhXMe9t5qPZLdkdRV6BldBscpfaxYJ0+kkG5jxb9YmjJEiAznJvb6vPzi7K7WrQF4ItF2nVTuY8W/WLof+vfpFSHd6hVy+4k6kp0a1oHOV2epTv1zlzlPlr0i6F1waOIuG6x9u8u4gIChNuTt3P4i/fIysp/eaVcoUW/mFm34wBZEdtpUlH75xcHfXtFcPgw/Pqr3UlUcaFFv5j5ZN4CAG5r0trmJMod6jXfTcBdt/Ph90vsjqKKCS36xcysbXOR02W5p10Tu6MoN7iqTASO2lOZt+0XHI78l1cqP1r0ixFjYF+yEJvRg5BgHYmyOIgKjaJy8DUcj1jO0qV2p1HFgRb9YmTDBjj99Tj+ff04u6MoN2pTMwGilzF1mt6kpa6cFv1iZObsswB06mRzEOVWLWMTIPwA38xKxmjdV1dIi34x8kbyLUTc15foaLuTKHdqFt2M2JB49h45wurVdqdRRZ0W/WLicPopjpb+hbhKl4w1r4q4plWbkjh0BYGHGzB1qt1pVFHnysDooSKyXETWishGEXneaq8uIstEJElEJolIiNVewnqfZM2PzbGtZ6z2rSLS2WPflR8aOfMXCMqgd0P9WIujcuWg7Y2GqVPRUzzqirhypJ8B3GSMaQA0BLqISDPgdeBdY0wtIBUYbC0/GEi12t+1lkNE6uEcY/daoAswUkS0i4mbfLtuDpwN5YEu2j+/OHr7t7dZ2TaarX842LTJ7jSqKHNlYHRjjDlhvQ22Xga4CZhitY8HelrTPaz3WPPbi4hY7V8bYzKMMTuBJHIZblEVzqbMOZQ70YayEWF2R1EeUDasLOmOfVDuD6ZNszuNKspcOqcvIoEisgY4BMwDtgNpxphzTwTZC5w7mVwV2ANgzU/HOXD6+fZc1sm5r6EikigiiSkpKQX+hvzR7t2GzEV/p0eVR+yOojwkoarz+Ciu3XI9r6+uiEtF3xiTbYxpCETjPDqv46lAxpjRxph4Y0x8hQo66pMr5s4VSHyQ4TffancU5SF1ytchPCScSo2Ws3YtbNdH7KtCKlDvHWNMGrAAaA5EiUiQNSsaSLamk4EYAGt+JHAkZ3su66gr8Pnin7kqLpl69exOojwlMCCQ+CrxHI9YAaBH+6rQXOm9U0FEoqzpMKAjsBln8b/dWmwg8L01Pd16jzX/Z2OMsdr7Wr17qgNxwHI3fR9+60xmFr9W7k3pHs/po5SLuYENBtKnfk/i47Xoq8ILyn8RKgPjrZ42AcA3xpgfRGQT8LWIvASsBsZay48FPheRJOAozh47GGM2isg3wCYgC3jYGJPt3m/H/0yYvwITmsbNVbWrZnE3qOEgAKQ3PPMM7NkDMTGXX0epi+Vb9I0x64BGubTvIJfeN8aYM8AdeWzrZeDlgsdUeflq2RxAGNalvd1RlBekn0mnZdcT8ExVvv0W/vY3uxOpokbvyC3iVqbNpVR6U+Kiy9kdRXlBnY/qMHbnv7j+ej3FowpHi34RtuvAMU5ELqdxhJ7a8RfxVeJZnryc225zjqZ18KDdiVRRo0W/CFuxOALe3cUTbYfZHUV5SUKVBLYc3kLn7scwBr77zu5EqqjRol+EzZkDkQFV6dbmKrujKC9JqJqAwXA6aiVxcXqKRxWcFv0iyuEwTDw5mOt6zCHIlT5YqliIrxIPwIp9y+ndGxYsgKNHbQ6lihQt+kXUzBVbOFl7HLENd9kdRXlRuZLlmNBzArfVvY3evSErC2bMsDuVKkq06BdRYxfOAeDBjjpMlr/p36A/ceXiaNIErr5aT/GogtGiX0QtPjCH4GPX0Oq6WLujKC87cuoIE9dP5OjpI9x2G8ydC8eP251KFRVa9IugtBNnOFzqF+oGa1dNf7T1yFbunnY3i3cvpndvyMiAH3+0O5UqKrToF0EzFu6DlHr0ur6L3VGUDRpd1YhACWR58nJatICrrtJTPMp1WvSLoHW/1CD400T+2bOb3VGUDcKCw6hfqT4r9q0gIAB69YKZM+HUKbuTqaJAi34RNGtuJq1aQalSdidRdmlapSkr9q3AYRz07u0s+HPn2p1KFQVa9IuYVdv2sfGWMlTu8I3dUZSNEqomkHYmje1Ht9OmDZQtq6d4lGv0tp4iZuScuRByiu4tr7E7irJR73q96VSzE9ER0YhAjx4wbRpkZkJIiN3plC/TI/0iZt7OOQScqkTvVvXtjqJsFBUaRUxkDGKNnHPHHZCerr14VP606BchWdkO9gTPIza7E0GB+k/n76ZsmsJzC54DoGNHqFoVPvnE5lDK57kyXGKMiCwQkU0islFE/m61lxWReSKyzfpaxmoXEflARJJEZJ2INM6xrYHW8ttEZGBe+1S5+2rBKkzYETrX0v75CpbuXcobS94gMzuToCAYPBhmz4bdu+1OpnyZK4eLWcA/jDH1gGbAwyJSD3gamG+MiQPmW+8BuuIc/zYOGAp8DM5fEsAI4AacI26NOPeLQrlm/dJK8PMLPNylo91RlA9IqJpARnYG6w+uB+C++5ztY8deZiXl9/It+saY/caYVdb0cZyDolcFegDjrcXGAz2t6R7ABOO0FIgSkcpAZ2CeMeaoMSYVmAfo3UUFsGxuDI2O/4drYyvaHUX5gISqztFKV+xbAUC1atCli7PoZ2XZmUz5sgKdGBaRWJzj5S4DKhlj9luzDgCVrOmqwJ4cq+212vJqv3gfQ0UkUUQSU1JSChKvWNt3+ARLUn6gXeeTdkdRPqJaZDXKlyzP8uTl59uGDoXkZJg1y8Zgyqe5XPRFJByYCjxmjDmWc54xxgDGHYGMMaONMfHGmPgKFSq4Y5PFwv/9OB9H31uJuWGF3VGUjxARmkc3Jz0j/XzbzTdD5cowerSNwZRPc6noi0gwzoL/pTFmmtV80Dptg/X1kNWeDMTkWD3aasurXblg+qbZkBnOkM4t7I6ifMh3fb9j6p1/3ZUVHOw8tz9zJuzZc5kVld9ypfeOAGOBzcaYd3LMmg6c64EzEPg+R/sAqxdPMyDdOg00B+gkImWsC7idrDaVj70px9gUMJGY090ID9M7b9RfAuTSH+HBg8EYGDfOhkDK57lypN8S6A/cJCJrrFc34DWgo4hsAzpY7wFmAjuAJOATYBiAMeYo8CKwwnq9YLWpfAwd/T8ITeelbv+0O4ryMRlZGXSY0IGRK0aeb6te3dlvf+xYyM62MZzySeI8He+b4uPjTWJiot0xbHX6NEQ90pnSpeHwe/qHkbpUjfdr0KRKEybfMfl829SpcPvtzjt0u+nDWP2OiKw0xsTnNk9v6/Rxn30GmeNmM+7miXZHUT4qoWoCK5IvvMDfvTtUqqQXdNWltOj7sDOZWbz2bjrNmgm3dihrdxzlo5pWacqu9F0cPHHwfFtwMNx7L/zwg7MLp1LnaNH3YY+P/Zrdva9mwPCtWM/VUuoSF9+kdc6QIc5z+p9+akcq5au06PuorGwHn257jRIZV3P/bXF2x1E+rHHlxnSs0ZHQoNAL2mvWhA4dYMwYcDhsCqd8jhZ9H/XclzPIiNzI4Gue1idqqssqFVKKuf3n0qFGh0vm3X8/7NoF8+bZEEz5JK0mPsjhMHyw+hWCjlfn7fv62B1HFREnM09ycW+8nj2hQgW9oKv+okXfB43+YTUno5ZzR5UnCQ3Rwc1U/iZvnEzEaxFsT91+QXtICAwaBNOnw/79ua+r/IsWfR805YPGlJu6gpFDB9kdRRUR11e6Hodx8NOOny6ZN2SI86mbn33m/VzK92jR9zHLlxvmz4en+scTFR6a/wpKAbXL1SY2KpbZSbMvmXfNNdCunXNULb2gq7To+5jbJvalRPfHefBBu5OookRE6FKzC/N3ziczO/OS+UOHws6dMH++DeGUT9Gi70OmL91EctQ3JDQsTenSdqdRRU3XuK6cyDzBkt1LLpnXqxeUK6cXdJUWfZ/y+JTXIbMknwz5m91RVBF0U/WbeKfTO9QuX/uSeSVKwMCB8N13cPDgpesq/6FF30cs3vAnO0p9SSPHUGrHlLc7jiqCwkPCebz541QpXSXX+fff77ygO358rrOVn9Ci7yOGffE2mABGDfqH3VFUEXYs4xhfrf/qgufwnFOnDrRpoxd0/Z0WfR9w8CD88ekT3Jj2GU1rR9sdRxVhu9J2cc+0e5jxx4xc5w8dCklJsHChd3Mp36FF3we89x5kplRj1CN32x1FFXHXVbyOqqWrMisp95HRe/eGMmXg73+H33/3cjjlE1wZLnGciBwSkQ052sqKyDwR2WZ9LWO1i4h8ICJJIrJORBrnWGegtfw2ERmY27780a6Daby15zY69lvPNdfYnUYVdSJC11pd+WnHT5zNPnvJ/NBQ5zn9I0egRQvo1w/27rUhqLKNK0f6nwFdLmp7GphvjIkD5lvvAboCcdZrKPAxOH9JACOAG4AEYMS5XxT+7v7RI8mK+5ZB92XZHUUVE13junIs4xi/7839UP7WW+GPP+DZZ2HKFKhdG154AU6d8nJQZYt8i74xZhFw8Vi2PYBzfQDGAz1ztE8wTkuBKBGpDHQG5hljjhpjUoF5XPqLxO8cSj3JTyfeo0JaV+66sZHdcVQx0b56e4ICgvh9T97nb8LD4eWXYfNm53CKI0Y4L/R+/bVzUHVVfBX2nH4lY8y5xzcdACpZ01WBPTmW22u15dV+CREZKiKJIpKYkpJSyHi+L/NsNg1f6IcJO8yLHf9tdxxVjESGRrLrsV081eqpfJetXh0mT3Ze2C1bFu66C1q3hpUrPZ9T2eOKH+FojDEi4rZjA2PMaGA0OAdGd9d2fYkx0PnZT9gf9R23hb3PA91a2B1JFTN59dXPS9u2zkI/bhz861/QtKnz6ZyvvOIca/f0aefpn5MnnV/PvXK+b97cOXCL8m2FLfoHRaSyMWa/dfrmkNWeDMTkWC7aaksGbryofWEh913kvfsuLHxnMN2eKMvUEXfaHUcVQ8cyjjF0xlB61+3NHdfe4dI6gYHOG7juvBNefBE++MB50dfVPv0REfDjj9Cq1RUEVx5X2KI/HRgIvGZ9/T5H+yMi8jXOi7bp1i+GOcArOS7edgKeKXzsouvfn87l5eca0LtXJb55VQu+8ozSIaVZtGsRgMtF/5zISHjrLXjgAWfRDwyEkiWdr1Klcp/OyoK774ZOneD776FjR098V8od8i36IjIR51F6eRHZi7MXzmvANyIyGNgFnKteM4FuQBJwCrgXwBhzVEReBM6N3PyCMebii8PF3sgflvDyju6U79ebz9/9kgC9S0J5iIjQpVYXvtvyHVmOLIICCn58FxcHL73k+vKLFkHnznDLLc4Lwr16FXiXygvk4uHVfEl8fLxJTEy0O4ZbzFqxlZuntiDobHk2Pv4bcdHl7I6kirlvNn5Dnyl9WHLfElrEeOe6UWqqszfQihXOQVv69fPKbtVFRGSlMSY+t3l6rOkFG3YepPukrmACmdN/lhZ85RUda3QkQAJyHVjFU8qUcQ7C3rYt9O8PH3/stV0rF2nR97CTJ6Htq4+TFXqAse1/oF3DGnZHUn6iTFgZ+tXvR8VSFb263/Bw5wXdW2+FYcPgjTe8unuVDx1124OysqBvX0hd8AEvjRrKvZ0S7I6k/Mz4nvY8Rzk0FKZOhQED4Kmn4NgxZ48gEVviqBy06HuIw2Ho/OQX/DyrDx++X56H77nR7kjKT53NPkvamTQqlKrg1f0GB8MXX/x19++xY86HC2oHBntp0feQW159k58jn6Ljk6d4+OEH7I6j/JQxhmtHXkvTqk358rYvvb7/wEDnEI0REfDOO3D8OIwZ42xX9tCi7wH3fjCOWVlPEZPeh5lv3m93HOXHRIRm0c2YuW0m2Y5sAgO8X21FnP3+IyLgv/+FEyfgq6+cfwko79M/tNxo458pRA+/nc9SBxOV2o51L3xGUKB+xMpeXWp14cjpI6zcb98DdUScD3V7+23nkz3794fsbNvi+DWtSG5gjPPcZas2mSQHLqZz4Kvsf2MuUeGhdkdTik41OyEIs7blPrCKNw0f7uzNM2mS845fHbbR+7ToX6HEP5Kp/eBz9B/goE7VqqwcsIPZ/36a0BA9c6Z8Q/mS5UmomsDs7d7rr385//wn/Oc/MHas85eAD98fWixpZSokh8Mw5KNP+XTfcKiQyROv38Frw68nMLCk3dGUusTLN71MSGCI3THOe/55Z2+e9993nut/4QW7E/kPLfqF8Pum3dw6eihHyswh8nQbpgwYS4fGteyOpVSe2tdob3eEC4g4nzZ74oSz/37p0s6/AJTnadEvAIcDRo0yPLLxFhyRO+gT/hFf/PtBvVirioRfd/3K3mN7uev6u+yOAjgL/6hRzsL/5JPOwv/gg3anKv606Lvg2MkMXps6m4WfN+P3nyrRtOcY3nmwIq2ui7U7mlIu+zjxY+bvnE+f6/oQIL5xoBIYCJ9/7nxcybBhzhu59CFtnuUb//I+6ExmFq9MmkvcE/cS+VIlXt3Zk3XZk/nkE1g2LUELvipyutTqwqGTh1hzYI3dUS4QHOwcsrFdO+doXd9+a3ei4k2P9HNwOGDxYvj865OMDa+JKXUQQiKolXkbgxr25fGnbqKk9sJURVTnmp0BmLVtFo0rN7Y5zYVCQ/8afKVPH5gxw/lsfuV+fv88fYfD8Pn8RD5Y8DVb/kzj1MSxhIVBjfue59amDXiqdxftb6+KjfjR8YQFh/Hrvb/aHSVXqalw002wdSvMmeMcpF0V3OWep+/1I30R6QK8DwQCY4wxr3lyf4fTT7Hijz2s37WHrfv3cvXRgezdKyzKeos/y4wjM2wPhJyAoGAqV7yVT75y0P3WAMLDR3gyllK26FKrC6NXjuZM1hlCg3zvYKZMGWexb9MGbr4ZXn3V+VcAOPvznztGPTed85g1JARKlPjrldf78HCIinJu1x+f+unVI30RCQT+ADoCe3EOn3iXMWZTbssX9kh/7Vq4+fmP2HfNCEzYkQtnvpFCpdLlCW05mrPV5lCxRAyNqjTgP7f3pHrlMrlvUKli4kTmCcKCwmx5Bk9B7N3rHIhlxw7P7SMkxDkecFTUha+cbRERl3+VLOmbvzh86Ug/AUgyxuwAsAZQ7wHkWvQLKzISqobWItLcSXRADDXKx1CnSgwNYmNIeLIMJcMAhlovpfxHeEg4AI/Nfox5O+ZdMK9yeGV+GvATAENnDGXJniUXzK9Vthbf9/0egH7T+rH6wOoL5jeo1ICven8FQK9JvfjjyB8XzG8e3Zwx3ccA0OnzTiQfT75gfofqHXi/6/sA9JnbkhLD04jL+mt+p2o9eDr+FUSg07SGZDnOQo6Ce2vsXQyt829OnjnL7fMa4jAX/kXQLnIIHcMf5/DxdF5NaYEjG7Ky4aAD9mdDxZ1/o8SGBziSuY+Urh0x+4B9OQIufgbW9YOy26Bvz/PNAYHOwh++7CVC/+xFVvk1pLe/58IPXiBq6buEJncio+JvpLa+Hy463o5YNIrg/a3IiJ7HySYvci+LGDUKt/N20a8K7Mnxfi9wQ84FROR8Nb766qsLtZPYWFj2VWdArwQplZuYiBjqVah3QVv5sPLnp6+OvJrUM6mXrHNObFQsGdkZF8yvHlX9/HSNqBqXDMZeLbLa+em4snFEhkZeMD86Ivr8dO1ytTmeefyC+XWqVCHaWqR+lbpkObIumF83ujLXXANZDqHx1gu/N4BOtStyT304nhHAkumXzu87sDy960HKyWCGzayHwwFZZ+FsFpw9Cx07lOXaENhxpASf7avH2Sxr/llwGKiXEEWVBpAWGMaKsHrni/q52n5do9JUqAdHQ0qxtoRz/+d/Zwk0ahFO+Sw4FBLBlpK1aRp3SUS38PbpnduBLsaYIdb7/sANxphHclu+OA2MrpRS3uJLA6MnAzE53kdbbUoppbzA20V/BRAnItVFJAToC0z3cgallPJbXj2nb4zJEpFHgDk4u2yOM8Zs9GYGpZTyZ17vp2+MmQnM9PZ+lVJK6bN3lFLKr2jRV0opP6JFXyml/IgWfaWU8iM+/ZRNEUkBdrmwaHngsIfjXAlfzqfZCs+X82m2wvPlfK5mq2aMqZDbDJ8u+q4SkcS87j7zBb6cT7MVni/n02yF58v53JFNT+8opZQf0aKvlFJ+pLgU/dF2B8iHL+fTbIXny/k0W+H5cr4rzlYszukrpZRyTXE50ldKKeUCLfpKKeVHimTRF5E7RGSjiDhEJM/uSyLSRUS2ikiSiDztxXxlRWSeiGyzvuY6+K6IZIvIGuvl0UdM5/dZiEgJEZlkzV8mIrGezFPAbINEJCXHZzXEi9nGicghEdmQx3wRkQ+s7OtEpLEPZbtRRNJzfG7PeTFbjIgsEJFN1s/q33NZxs7PzpV8tnx+IhIqIstFZK2V7flclin8z6sxpsi9gLpAbWAhEJ/HMoHAdqAGEAKsBep5Kd8bwNPW9NPA63ksd8JLefL9LIBhwP+s6b7AJB/KNgj40Kb/a22AxsCGPOZ3A2bhHPmuGbDMh7LdCPxg0+dWGWhsTZcG/sjl39XOz86VfLZ8ftbnEW5NBwPLgGYXLVPon9cieaRvjNlsjNmaz2LnB2E3xmQC5wZh94YewHhrejzQ00v7zYsrn0XOzFOA9iIieJ6d/075MsYsAo5eZpEewATjtBSIEpHKPpLNNsaY/caYVdb0cWAzzjGyc7Lzs3Mlny2sz+OE9TbYel3c46bQP69Fsui7KLdB2L31j1rJGLPfmj4AVMpjuVARSRSRpSLS04N5XPkszi9jjMkC0oFyHsxUkGwAva1TAFNEJCaX+Xax8/+ZK5pbpwlmici1dgSwTj00wnnEmpNPfHaXyQc2fX4iEigia4BDwDxjTJ6fXUF/Xr0+iIqrROQn4KpcZv3LGPO9t/Nc7HL5cr4xxhgRyatfbDVjTLKI1AB+FpH1xpjt7s5aDMwAJhpjMkTkAZxHODfZnKkoWIXz/9gJEekGfAfEeTOAiIQDU4HHjDHHvLlvV+STz7bPzxiTDTQUkSjgWxG5zhiT67WbgvLZom+M6XCFm/DoIOyXyyciB0WksjFmv/Xn6qE8tpFsfd0hIgtxHm14oui78lmcW2aviAQBkcARD2QpcDZjTM4cY3BeM/EVHv1/diVyFjFjzEwRGSki5Y0xXnmYmIgE4yyoXxpjpuWyiK2fXX757P78rP2micgCoAuQs+gX+ue1OJ/esXMQ9unAQGt6IHDJXyYiUkZESljT5YGWwCYP5XHls8iZ+XbgZ2NdJfKwfLNddJ63O87zr75iOjDA6onSDEjPcWrPViJy1bnzvCKSgPPn3Ru/yLH2OxbYbIx5J4/FbPvsXMln1+cnIhWsI3xEJAzoCGy5aLHC/7x6+8q0O15AL5zn/zKAg8Acq70KMDPHct1wXpXfjvO0kLfylQPmA9uAn4CyVns8MMaabgGsx9lbZT0w2MOZLvksgBeA7tZ0KDAZSAKWAzW8+Hnll+1VYKP1WS0A6ngx20RgP3DW+j83GHgQeNCaL8BHVvb15NGbzKZsj+T43JYCLbyYrRXOi4/rgDXWq5sPfXau5LPl8wPqA6utbBuA56x2t/y86mMYlFLKjxTn0ztKKaUuokVfKaX8iBZ9pZTyI1r0lVLKj2jRV0opP6JFXyml/IgWfaWU8iP/D0C1CFQcw0pLAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Within snapshot region, find stars within 30 kpc:\n", "ind_in_snap_30kpc, in_snap_30kpc = hy.crossref.find_id_indices(ids_stars_30kpc, stars.ParticleIDs)\n", "\n", "# Verify that we found all of them\n", "print(f\"Located {len(in_snap_30kpc)} out of {len(ids_stars_30kpc)} particles within 30 kpc in snapshot region.\")\n", "\n", "# Make a radial histogram\n", "radii_30kpc = np.linalg.norm(stars.Coordinates[ind_in_snap_30kpc, :] - centre[None, :], axis=1) * 1e3 # Convert radii to kpc\n", "hist_30kpc, edges = np.histogram(np.log10(radii_30kpc), bins=30, range=[-1, 3])\n", "\n", "ax = plt.gca()\n", "ax.plot(midpoints, hist, color='blue')\n", "ax.plot(midpoints, hist_30kpc, color='green', linestyle='--')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, we now (green) get zero stars beyond 30 kpc (the truncation does not coincide with a histogram bin edge, so it appears smooth here...).\n", "\n", "Note that another approach to get the same result here would have been to only load a snapshot region extending to 30 kpc from the centre, i.e. with" ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Checking 1 cells...\n", "Reading 'Coordinates' took 0.112 sec.\n", "Region setup took 0.210 sec.\n", "Selection region contains 1 cells, 3 segments, 859127 particles, 3 files\n", "Exact selection region contains 80110 particles.\n" ] } ], "source": [ "stars_30kpc = hy.ReadRegion(snap_file, 4, centre, 0.03, exact=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Off-centre FOF haloes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Working with FOF haloes that are not the central cluster is quite straightforward in principle -- just remove the condition of targeting (satellites in) the first, or most massive, FOF group:" ] }, { "cell_type": "code", "execution_count": 95, "metadata": {}, "outputs": [], "source": [ "# Read M200, Mstar, and central galaxy ID for all galaxies in snapshot 29\n", "m200_snap29 = hy.hdf5.read_data(sim.fgt_loc, 'M200', read_index=29, index_dim=1)\n", "central_galaxy_snap29 = hy.hdf5.read_data(sim.fgt_loc, 'CenGal', read_index=29, index_dim=1)\n", "mstar_snap29 = hy.hdf5.read_data(sim.fgt_loc, 'Mstar', read_index=29, index_dim=1)\n", "ngal = len(m200_snap29)" ] }, { "cell_type": "code", "execution_count": 101, "metadata": {}, "outputs": [], "source": [ "ind_satellites = np.nonzero((mstar_snap29 > 9.0) & (central_galaxy_snap29 != np.arange(ngal)) &\n", " (m200_snap29 > 13.0))[0]\n", "ind_central_satellites = np.nonzero((mstar_snap29 > 9.0) & (central_galaxy_snap29 == 760) &\n", " (central_galaxy_snap29 != np.arange(ngal)))[0]" ] }, { "cell_type": "code", "execution_count": 102, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Found 193 satellites above 10^9 M_Sun in total, 69 in central cluster.\n" ] } ], "source": [ "print(f\"Found {len(ind_satellites)} satellites above 10^9 M_Sun in total, \"\n", " f\"{len(ind_central_satellites)} in central cluster.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, there is a caveat! Not all these outer groups and their satellites are necessarily (well) within the resolved region of the simulation. To filter out bad ones, look at the `ContFlag` data set. This is 0 for any galaxy >= 8 cMpc from the nearest boundary particle, 1 (2, 3) for >= 5 (2, 1) Mpc distance, and 4 if it is < 1 cMpc away." ] }, { "cell_type": "code", "execution_count": 104, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Out of 193 satellites, 125 are far from the boundary.\n", "53 are dubious.\n", "15 are very close to the boundary.\n" ] } ], "source": [ "contflag_snap29 = hy.hdf5.read_data(sim.fgt_loc, 'ContFlag', read_index=29, index_dim=1)\n", "contflag_all_sats = contflag_snap29[ind_satellites]\n", "print(f\"Out of {len(ind_satellites)} satellites, {np.count_nonzero(contflag_all_sats <= 1)} are far from the boundary.\")\n", "print(f\"{np.count_nonzero((contflag_all_sats > 1) & (contflag_all_sats < 4))} are dubious.\")\n", "print(f\"{np.count_nonzero((contflag_all_sats >= 4))} are very close to the boundary.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The central cluster's satellites should, by definition, all be very far from the boundary:" ] }, { "cell_type": "code", "execution_count": 106, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 central cluster satellites are not far from the boundary.\n" ] } ], "source": [ "print(f\"{np.count_nonzero(contflag_snap29[ind_central_satellites] > 0)} central cluster satellites are not far from the boundary.\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.6.8" } }, "nbformat": 4, "nbformat_minor": 2 }