Skip to content

Instantly share code, notes, and snippets.

@aphearin
Created April 7, 2018 04:37
Show Gist options
  • Save aphearin/9052a6192ff925e5e430533f8e583e97 to your computer and use it in GitHub Desktop.
Save aphearin/9052a6192ff925e5e430533f8e583e97 to your computer and use it in GitHub Desktop.
Clustering and lensing of isolated vs. non-isolated centrals
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 153,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"morange = u'#ff7f0e'\n",
"mblue = u'#1f77b4'\n",
"mgreen = u'#2ca02c'\n",
"mred = u'#d62728'\n",
"mpurple = u'#9467bd'"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1444262\n"
]
}
],
"source": [
"from halotools.sim_manager import CachedHaloCatalog\n",
"halocat = CachedHaloCatalog(simname='bolplanck')\n",
"print(len(halocat.halo_table))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"from halotools.empirical_models import PrebuiltSubhaloModelFactory\n",
"model = PrebuiltSubhaloModelFactory('behroozi10', redshift=0.)"
]
},
{
"cell_type": "code",
"execution_count": 91,
"metadata": {},
"outputs": [],
"source": [
"model.populate_mock(halocat)\n",
"\n",
"from halotools.utils import crossmatch\n",
"idxA, idxB = crossmatch(\n",
" model.mock.galaxy_table['halo_id'], halocat.halo_table['halo_id'])\n",
"\n",
"model.mock.galaxy_table['halo_halfmass_scale_factor'] = 0.\n",
"model.mock.galaxy_table['halo_delta_vmax_behroozi17'] = 0.\n",
"model.mock.galaxy_table['halo_nfw_conc'] = 0.\n",
"\n",
"model.mock.galaxy_table['halo_halfmass_scale_factor'][idxA] = (\n",
" halocat.halo_table['halo_halfmass_scale_factor'][idxB])\n",
"model.mock.galaxy_table['halo_delta_vmax_behroozi17'][idxA] = (\n",
" halocat.halo_table['halo_delta_vmax_behroozi17'][idxB])\n",
"model.mock.galaxy_table['halo_nfw_conc'][idxA] = (\n",
" halocat.halo_table['halo_nfw_conc'][idxB])\n",
"\n",
"mask = model.mock.galaxy_table['stellar_mass'] > 10**9.5\n",
"mask *= (model.mock.galaxy_table['halo_mvir_host_halo'] > 0.)\n",
"galaxies = model.mock.galaxy_table[mask]"
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['halo_upid', 'halo_hostid', 'halo_mpeak', 'halo_x', 'halo_y', 'halo_id', 'halo_z', 'halo_vx', 'halo_vy', 'halo_vz', 'halo_rvir', 'halo_mvir', 'halo_mvir_host_halo', 'x', 'y', 'z', 'vx', 'vy', 'vz', 'galid', 'stellar_mass', 'halo_halfmass_scale_factor', 'halo_delta_vmax_behroozi17', 'halo_nfw_conc']\n"
]
}
],
"source": [
"print(galaxies.keys())"
]
},
{
"cell_type": "code",
"execution_count": 123,
"metadata": {},
"outputs": [],
"source": [
"from halotools.mock_observables import conditional_spherical_isolation\n",
"\n",
"mstar_mask = galaxies['stellar_mass'] > 0\n",
"x = galaxies['x'][mstar_mask]\n",
"y = galaxies['y'][mstar_mask]\n",
"z = galaxies['z'][mstar_mask]\n",
"sample1 = np.vstack((x,y,z)).T\n",
"sample2 = np.copy(sample1)\n",
"marks1 = galaxies['stellar_mass'][mstar_mask]\n",
"marks2 = galaxies['stellar_mass'][mstar_mask]\n",
"\n",
"r_max = 1.5\n",
"cond_func = 2\n",
"galaxies['is_isolated'] = conditional_spherical_isolation(\n",
" sample1, sample2, r_max, marks1, marks2, cond_func, period=model.mock.Lbox)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 124,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"isolation fraction = 0.40\n"
]
}
],
"source": [
"mstar_mask = (galaxies['stellar_mass'] > 10**10) & (galaxies['stellar_mass'] < 10**10.25)\n",
"print(\"isolation fraction = {0:.2f}\".format(galaxies['is_isolated'][mstar_mask].mean()))\n",
"sm_sample = galaxies[mstar_mask]\n",
"\n",
"cenmask = sm_sample['halo_upid'] == -1\n",
"isomask = sm_sample['is_isolated'] == True\n",
"all_iso = sm_sample[isomask]\n",
"all_cens = sm_sample[cenmask]\n",
"all_sats = sm_sample[~cenmask]\n",
"iso_cens = sm_sample[cenmask & isomask]\n",
"iso_sats = sm_sample[~cenmask & isomask]\n",
"noniso_cens = sm_sample[cenmask & ~isomask]\n",
"noniso_sats = sm_sample[~cenmask & ~isomask]\n"
]
},
{
"cell_type": "code",
"execution_count": 209,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEhCAYAAABr1YsqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHBlJREFUeJzt3bGPImu+3vHnZ21kJxx2N5nE5xaWHG3CdEsOVt5gqhNLJ7kLZ7KbXA1YsjZtdBzcHd3AIyaydH0lw9H9A+bAbnJkJ7CWbG1wpekmvNJaatYOth2sh6nEwSb3dcBbTDVTQNEU0Lx8PxKagSqol+ruh5dfvVWvOecEAAjHPzl2AwAA5SLYASAwBDueFDOrm1nXzK6P3RZ8YmaRmcXHbgeKIdjPmA/RwYplLTNr+NvBQtY5N5H0XtLVIbbnP0TuzMyZWWXNepGZffTrdteteyhl//z8c1pm1vO37HusSxr4/fTRzEZmVi/jfaB8Pzh2A3B4/g/ypb8b5SxvSZJzbujvR2bWc861D9TEyYG2I+dcx8w+SOpqvi9WbTuWVJH0wn/4HM0+fn5m1nLO9Zde41ZSLX3MOfeFmVWcc0k57wT7Qo/9DDnnJs65jqR3K1ZpZ//InXNTzYNtJ743ONr1dfZg6m+fhaQkmVlD0kxafKPYi6L7p+yfX963D//86nL5hVA/DQQ7HvB/5HlfsZPH1FjNrGJm1z6wZs65g5RYivLvaaJ5sF/mLI8yy8Z72H6p++eRP79I0nLpRVrzYYenjVIMlkWS8nplM80Do1C4+ZD4xr9ezzn3dpdG+dLALNPGftp7zGzrvV829etdFSgfRc65sZmtCrG6c27oa9m9Xd5DVtn7J2Prn59zbmJmz3N649l9qaUPhroyPwM8LQQ7llX1KUCzEkk/3PRk38Pt+Nd5U0bpwsx6krq+pJCG4kCfDrB+q3k4jv3yO+dcTdJwi83cSbpY2m4saey3F6mEHvs+9s+SR/38ltvhy0/TdJ/KH3vI/AymevgzwBNCsKMUPrC6mofKIoRLeN26pIvs6znnEjObZg74xZJeLT9vU2j6wE5D8EGPPS1L+G01/P8fHcL72j/7kPk28SJ9bLm9zrmpPyi7cT/j8KixI08157GKpA8Fnpsov8f4WBfKlAMy7iQ99/+f6WGbq8ovRyyL9akXPtH8PS6WZXqrVyqvvl72/smzy89Pmn8ANQuUWRItfcvB00CwY9mNHgZcqqo1wxCdc1PnXFPzOnTXj4Pe94G3NMB68qM+/DZvCvaIq2l4ZUoM6Yk42SCPJe00mueA++dRP7+UH/P+4BuF3yd5F5Waaf8fUngEgh0P+KCb5oyQqGR6sOueP/UHLDuS2mY22PGMxRvlj/Ko6VPYTiXNfMkk3mFkSZJua+nAbG593QdeY5sN7GH/LL/+o39+/gD1cCnUY83DO+8g9IUOeM4BiiPYz1veV3Zp/lX8m/SOr3NvVYpwziV+rPUrSXV/puJWIehfZyJpkj3L0YfWRWas9qVzbuhv/dwXWrLirMmppJdLARhLSnIOLl5n/t/Ytvdd0v551M/PfyANsuHvA/wm882lkn7g5JVk/IfAd0/5WMFZc85xO7Ob5j3QruZnFjrNywOtpXVamodaQ9J1SdttSRoVaNvAt+s68/i1b0vD/7+SWRZL+ujfz61/fmPNNnp+/bvsen6fRJl2pOul+yi7zbpff6D5t4SD7J8yfn6Z/ZV9r27FLfuer/3rXpf1O8FtPzfzPzDgJPmeckN+TLXvhVY1L3WkveJ9bLer+cHIdAjhG8eYbjwRBDtOmg/YnsspCZjZyO3pTNdM6SWW9J3E6fZ4Oqix49S9V851UHxdeW8H9vwHyUzzunRCqOMpoceOk+cP8mVPf69oPpSx0IFUIDQEOwAEhlIMAASGYAeAwBDsABAYgh0AAnP0y/b+6Ec/cl9++eWxmwEAJ+X29vb/Oud+nLfs6MH+5Zdf6ubm5tjNAICTYmb/e9UySjEAEBiCHQACQ7ADQGAK1dgz156+lPTebZhR3a8/lb9eNKd2A8DhbAx2M+u5+Ywv6f1bM9OqcPdX23vvnBum982skd4HAOzX2lKMv7b18lXresrMzpKjtRTiI+VPqwUA2INNNfaqpOucab/yJstdNd3YTDmXVQUA7MfaYPfXnH6+NInBlVbPf1nV57OWZycFBgDs2cZRMS4zia8P51irSyvptGRZadDnTrx7f38vM1vcXr9+vbHRAIDVtj3zdCDpRd40ZF7eLDJpoC/35CVJz5490/39/ZbNwLa++pvf5j7+/S9+euCWANi3wuPY/WiXbrYHn2Omz+vvFYn5IAHgUAoFu5k1JI2cc2N/P+8gaVq2WQ7wqlbX5IGgjcdj1Wo1tdvt3PvHVKRtzWZT/T6noZyaIuPYY/lw9jX2qqSX8hMF+xEz9cwQx/7SuPUrzYdIAlv5/c8bR9nun/2qvFMu4jhWp9PR7e1t7v2y9ft9tVqt0trWbrcVRQ8HxW2zDRzH2mD3QT7yd7PhnP3NjyU108eccx0zu/a9/EjSHScnAfuXJInu7u5Kfc04fjhSeR/bQPnWBruvi9uGdfqS+kuPrb3kAIByJUmiV69efda73vU1p9OpKpWKoijayzawH1wEDCjBZDLRcDjUcDhUu91Wkuw2VuDt27caj8eL19u0nfF4rCRJNJlMFs9NTadTdTodDYdDdTqdwm2bzWbqdDrqdrs7bWPVe8H+HH2iDSAEzWZTg8FA9Xpds9lMb968WQTitvr9vur1+qIMMp1+Gl28ajuNRkPT6VQfPnzQ9fX1g9e7urrS7e3touf96tUrDQaDje2IokjNZnNRc3/MNta9F+wPPXagBKPRSPX6fLDYxcWFJpN1o4LXi6JI7XZb/X5fSZI8OFC57Xb6/b6iKFKlMh+FXK/XH/S0y7BuG+veC/aHHjtQgiiKNBwONZvNlCSJZrPc8/EKieNY3W5XvV5P7XZbrVZLvV7vUdtJD3Rmw/zrr79+dNu23ca694L9occOlOD58+eKokitVuuzkSTbGo/HajQaGo1G+vjxo6bT6aKEUXQ7w+F8INrl5aUqlYriOF7cygrWIttY916wPwQ7sKP0oGJaIkmDKz3QuK3RaLR4XqVSWbzupu2kI1ey0rp49vFdSjHbbmPVe8F+UYoBdhTHser1+qLWnN76/b7iONZgMNB0OtVwOFQURQ/uNxqfn4RVq9Ue9GxrtdriNVdt5/r6WvV6Xe/evVssTw0GA71580aXl5eStFg2mUzWti2KIvV6PSVJovF4rDiO1Wg0ttrGqveC/TLn3FEbcHFx4W5ubo7ahnPARcCAsJjZrXPuIm8ZpRgACAzBDgCBIdgBIDAEOwAEhlExZ46DqkB46LEDQGAIdgAIDMEOAIEh2AEgMAQ7AASGUTF4uno/O8522//9ONs9gGazqaurq0dfF308Hqvdbpd6lcintL2yHav99NiBM5KGzGPFcaxOp7PVc/r9/uaVStxeWXZpd+pY7afHDpyRXa8Vv60kSRYTcZySU213ih47cCbS67YfaqKLJEn06tWrg2yrTKfa7ix67MCOxuOxOp2O4jheXI98NBqp0+ksrj2eJMniGuaz2UxRFCmO40LPLctsNlu8blrvffv2rer1upIk0Wg0Wjy+qr15sh8Wo9FI3W5XlUplMTHIZDJZbCc7qXWv19Pl5aXev3+vb775ZjFnapIkD67tXqTnnH1OpVJRtVp9MCFJ3rZW7fuf/OQnue1O6+XtdluVSkW9Xk+/+c1vVKlUVu6DPKv2eZkIdmBHcRyr3W6r1+up2+1Kmodo9v6LFy90e3u7eE6z2VS1Wi303LJEUaRms7loR7/f/yxsU6vamzcDUrPZ1GAwUL1e12w205s3b9TtdhczK3348EHX19cPnnN1daXb21tVKhVFUaRXr15pMBgstp0GpiS9f/9+43t78eKFBoOBoijSZDJRs9lcfCCs2taqfX93d6erq6vP2p2u/+7dO93e3qparW7cB8vW7fMyUYoBSpLtYVer1cVUcelsRFkvX77UmzdvNj53n6IoUrvdVr/fV5Iki5EyRdqbNRqNFoF/cXGxcTrA9JtAGtz1en0xlV76b7a3W6vV1r5e+py0zfV6XaPRaOO2Utvu+3T9RqOxeN2i+2DVPi8bPfZz8cff5T/+43952HYELNuDk+a9P2ne41xeln593/TcPP1+/0FvOk+RUk4cx+p2u+r1emq322q1Wur1eoXamxVFkYbDoWazmZIkWdt26VNpJRuwX3/9taR5WWd525vkPSd97+u2lcrb9+s+TPL2a9F9sGqfl41gB/asVqstepCp7KTU2yqrlzcej9VoNNRoNJQkiZrNpqbT6dbtff78ub799ls1Gg1NJhO9e/cud710jtfLy0tNp9MHNfv0/+m8rdtY95x12ypq1dy0WUX3wap9XvbxFEoxwJ61Wq3Paqnv3r3TN998c6QWzY1Go0UvvFKpLIJ7m/amB0izByqlTyNwoij6rLSR1t6zj6c96jR0s8s2fTuJ41hJkjxoczoGfd221slr9yqb9kHWqn1eNnrseLpO5AzQyWSiXq+nJEkW9eler6ebm5tFb28wGKjT6Sx6kO12W/V6vdBz99HO8XisWq2m6XS6CKJarbboOa5r72Aw0HQ6XbSvXq8vatnprd/v6/r6etGbTpenBoPBg5Evecuurq4W4frdd9/p+fPnK7+t3N7eqtPp6OrqStLDXvmqba3b9y9fvtRsNnvQ7vF4rHfv3ilJEtVqtUVb4jheuQ/iOH6wv9bt8zKZc670F93GxcWFu7m5OWobzsFXf/V3+QtW1NiZaAN42szs1jl3kbeMUgwABIZgB4DAUGMPzKo5TAGcD4L93K0a3y5q7MCpohQDAIEh2AEgMAQ7AASGYAdK0Gw2d5pxJz1hqN1ul9iqp7O9sp16+/eNYAdKwJRzxZ3ylHOnotCoGDNrSLp0zq3dk369SNJQ0kxSS9LQOXeYKVuwZpTL6TnW0M3HnHXLlHPFnGq7T83aHruZxWZ2LaktKX86kIeqkrqS7iT9XtKUUEfomHKumFNt9yla22N3zo0ljc3shyoW7JL0haQqgY5zwZRz5zPl3Kko/QQl51wiaf/TvwBPBFPOnc+Uc6ei9IOnZtYys4b/93rzM4CwMOWcHjxPCmfKuVNRdo99LGnme+0ys56ZtZxzux8GB04EU87NhTjl3KkotcfunJumoe6NJK0dSXN/fy8zW9xev35dZpOAg0unPxuNRvr48eNiYoVarfZZMG2aci6KIrVarbWjbobDoSQt6uBxHC9uabilpYxtrHvOum0VlbZ7naL7YNU+P1elBbuZVczMmVn2yEai+fDHlZ49eybn3OJGsOPUMeXcaqc+5dypKLsU83apxx5JOt+PTZwFppw7nynnTkWhqfHMrCup4pxrLz0eSao754b+/rVz7m1m+UhSL12eh6nxyrVyCrwtff/Xf1nK6+zilE5QAg5t3dR4a3vsZlaXFEtqSKqa2Z2ksXMu/R4US2pqfqapJPX9SJhEUk0bQh1Yh4AFHmfTCUoTSRNJb1cs70vqZ+4nq9YFABwGMygh17oyCD1p4Gnj6o4AEBiCHQACQ7ADQGAIdgAIDMEOAIEh2AEgMAQ7AASGYAeAwBDsABAYgh0AAkOwA0BgCHYACAzBDgCBIdgBIDAEOwAEhmAHgMAQ7AAQGGZQQr4//m7NQmZQAp4yeuwAEBiCHQACQ7ADQGAIdgAIDMEOAIEh2AEgMAQ7AASGYAeAwBDsABAYgh0AAkOwA0BgCHYACAzBDgCBIdgBIDAEOwAEhmAHgMAQ7AAQGIIdAAJDsANAYAh2AAgMwQ4AgflBkZXMrCHp0jnXKbDutaSppKokOef6O7UQALCVtT12M4t9ULclVTa9mJl1JU2dc0Mf6DX/oQAAOJC1we6cGzvn3kqaFHy9lnNumLk/0vxDAQBwIKXV2M2snvPwTFJc1jYAAJuVefC0qnmQZyWSZGYbyzgAgHIUOnhaUEX+gGlGGvRV+ZBHCXo/W7PwLw7WDABPU5nBnhfcadAv9+QX7u/vZWaL+7/85S/1+vXrEpuFsn31N7/Nffz7X/z0wC0BkKfMYJ/p85EzFUlyzq3srT979kz39/clNgMAzltpNXbn3ESf99qrksZlbQMAsNlOwW5m0dI49f7S/StJvV22AQDYztpSjB/CGEtqSKqa2Z2kse+dyy9rShpKknOuY2bXPtwjSXdL49oBAHu2Nth9gE8kvV2xvC+pv/RY7roIyB9/t2IBB0+Bp6DMg6c4kK/+wJBGAKtxdUcACAzBDgCBIdgBIDAEOwAEhmAHgMAQ7AAQGIIdAAJDsANAYAh2AAgMwQ4AgSHYASAwBDsABIZgB4DAEOwAEBiCHQACQ7ADQGAIdgAIDMEOAIEh2AEgMAQ7AASGYAeAwBDsABAYgh0AAkOwA0BgCHYACAzBDgCB+cGxG4BwfPU3v819/Ptf/PTALQHOGz12AAgMwQ4AgSHYASAw1NhRnj/+bsUCauzAIdFjB4DAEOwAEBiCHQACQ7ADQGAIdgAIDMEOAIEh2AEgMAQ7AASGYAeAwBQ689TMriVNJVUlyTnXX7NuQ1IkaShpJqklaeicm+7cWgDARht77GbWlTR1zg19oNd8eK9SldSVdCfp9/65hDoAHEiRUkzLOTfM3B9Jam94zheSas65L5aeCwDYs7WlGDOr5zw8kxSve55zLpGU7NAuAMAjbaqxVzUP8qxEksys4gP8M2bW8s+rSqo4597u2lAAQDGbSjEV+QOmGWnQLz+eGkv6bqkm31q1gfv7e5nZ4vb69esi7QYArLCpx57XI08DfbknL0nKOVA60vxgau5ImmfPnun+/n5DMwAARW3qsc8077VnVaRFHf0BM6uYmTOz7HMSzYc/AgAOYG2P3Tk3MbPlAK9qXm5Z5e1S6Eeaj4HHtno/W7HgLw7aDACnpcgJSn0za2SGLV5J6qULzSySVPc19cTMPiw9vympU05zgdP3+5/nnwbyZ78qb2Twqm2sU+b2cVwbg9051zGz68wZpXdLY9NjzcM7fazvz1RNJNUk9RjLDgCHU+iSAuuGK/qRL/3M/UQSwxuBA3hMzxzhKxTswC6++qu/y338+7/+ywO3ZLN1QbmqVHGI0gqwDYIdZ+kp9nQf86EC5CHYEbQyA3zb1zrmtnHeCHYEgeDbH0pNp4eJNgAgMPTYcVLome8P+zYc9NgBIDD02AGUipr88RHseHIoCZwGfk5PF6UYAAgMPXYAB8EJWIdDsONo+CoP7AfB/oR99Qeuu47zxoHYxyHYARwd397KxcFTAAgMPXYczbPa3+c+fn/3rw7cEiAs9NgBIDD02LF3//inP+U+/uf273Mf/0/6b/tsDhA8euwAEBiCHQACQykGQDAY9z5HsAM4ayF+GBDsT0HvZysWcOYpgO0R7CjNqtEvAA6LYAdwcrgEwXqMigGAwNBjx5Oz6lIDEpcbwOOcWw+fHjsABIZgB4DAUIrB1hj9AjxtBDuenFUXB5O4QBhQBMEOADlO+YxUgv2QVp5hCgDlIdhxUph1CdiMUTEAEBiCHQACQykGuRjSCJwugn0fOEi6N8yTCmxGsCMIHFQFPqHGDgCBKdRjN7NrSVNJVUlyzvXLXP8klVhu+eoPzJS0K0o0OJRTOHFpY4/dzLqSps65oQ/ompmtvAbmtuuH4PX3/+vYTXi0f/zTn3Jvu/ifv/0vJbVud89qf597O5T/+A//cLBtnQL2x0OvX7/ey+uac279CmYfnXNfZO7HkjrOuasy1r+4uHA3NzePavxBFOiZ27/9H3L/+V8/ehPH7LHvY/TLf3377/Rvrv+29Nct06/df8h9vOyafPTrX2n65z8v9TVPWcj74zE9djPTpgxe89xb59xF3rK1pRgzq+c8PJMUl7H+kxL4SBaGLz60qnTz61p+4EsciMXp2FRjr2oezFmJJJlZxTmX7Lj+bp5oGG/bAyd0n451V5bUv8h/+FC9f6CotaUYXxv/dqm0UpH0UVLNOTfdZX2//P9J+qeZh/6PpPvHvZ2jeabTa/M+sT8+YV88xP54aJf98c+dcz/OW7Cpx57Xw676f5d75o9ZX865f7ahDQCALWwaFTOTVFl6rCJJK8oq264PACjZ2mB3zk30eS+8KmlcxvoAgPIVOfO0vzQO/UpSL71jZtHS8rXrnzIza/hx+o9aHpp179fMrv1t4E9YC96q/WFmFb8vWmbWPYffkaJ/C2YWRDZssuZ3o+F/N6LM70m06/Y2nnnqnOv4jTUkRZLunHPZAZuxpKakYcH1T44fi1/X/EMq7wDw2uWhKbA/es65dub+rR+v+/aAzTyYAj//b5xzncz6t2bWCvGM7G3+FnzQ5Y7DDkWB/VGV1PW3RNKrvEEm2yp0SYF1f5D+l7O/9FhQf8DOubGksZn9UJ8fQ9i4PDTr3q8fBbVcjutp/osb1O9FqsDPv2Fmd5kgn2r+hx5csBf9WyijV3oKCu6PLyRVywj0FBcBQ9mqkvK+Tgb/gbfG1VLvPJL0/liNeSJiSaNjN+IpcM4lZYa6xGV7UTLn3NTMni/9ol7pjA+gZ/dFenZ2aN9qt+HLE98p8DJMUWbW0nxEYVVSpYzfDYIdpfOjoyQtSjOxpOfHa9Hx+f3wtebHo14duTnHVnHOJWZ27HY8BWNJs3Q4uJn1yjj+QikG+zaQ9KLsr5qnxn/d7vuL4X3re2lnx8wapz6YokzOuenSOT4jSZ1V6xdFsGNv/KiHbrYHf458bz2rp0CGAG/DH3fhREXPD290S78fiebHYHZCKQZ74Ye7jvyoAJlZ/RwD3teTR2b2xfLZ13u5MN7TVpcUZa4Ceymp4s9zGJ7pt7q3S78DkUoYMk2wo3Q+zKqaD/Oq+P+/lHR2wS7pRlJ/6Y/3SvMgO6dQ13IJxpejonM9kOyPM3xYeripEkoxBHsBvocRS2pIqprZnaRx2gPdtDw0696vD/J0GFu23BBsXXXd/vB/vL3M2bc/1HyGsZ3/eJ+ion8LPtSbmvfgr/X5h18QCuyPvn//iaSapF4ZxyA2zqAEADgtHDwFgMAQ7AAQGIIdAAJDsANAYAh2BMlf9/zOj0LYtJ7z/57FFQcRPoIdQfLDCTvSpwtvrTF1znXO9AQZBIhgR+iGktp5C9KzQg/bHGD/CHaErqf5VRWBs0GwI2i+vDJdmoc3vTAXpRcEiWDHOejp83LMRV5N3czqfs7ewjM+PeY5wD4R7Aien7QgLhK8/hoeLzW/cFnR159ofmGvws8B9olgx7kYSmpJi1EyN2vWfUyJJrgLWOF0Eew4F9lyTDXEKwkCKS7bi7PgnBubWZEx7anIn7B0pfmlVKfSYojkTPNyTW9FnT49UJtoPr9nsJcsxtNEjx0hWz6TdCjp23RWpw0qfr33ml9LO5VO9fdGOdPb+Q+OS+fc0D8/8h8GwMEQ7AhSOt+qn+QiPWja03xW+HSdhuZnp0b+kgLZg6vpRAiJ5pNjpF7456WzRC17qfmHQfZ1mju9GWBLlGIQpOwlBTKPTbOP+RJJ4TKJD/6Bc+7K3+8WfCqjZXBQ9NiB4r6W78n7kE9yyizvNK/Lp+r+MeBgCHYgw9fI65LaPrybmo+BjyR9J6niSzEXmpd1ouxzfP391swafr2Eg6c4NOY8BYDA0GMHgMAQ7AAQGIIdAAJDsANAYAh2AAgMwQ4AgSHYASAwBDsABIZgB4DAEOwAEJj/Dwkz039wpMEbAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x12354db90>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1, 1)\n",
"\n",
"logm_bins = np.linspace(10.5, 15, 50)\n",
"\n",
"__=ax.hist(np.log10(all_sats['halo_mvir_host_halo']), bins=logm_bins, normed=True, \n",
" alpha=0.8, label=r'${\\rm all\\ satellites}$', color=mred)\n",
"__=ax.hist(np.log10(noniso_cens['halo_mvir_host_halo']), bins=logm_bins, normed=True, \n",
" alpha=0.8, label=r'${\\rm non-isolated\\ centrals}$', color=morange)\n",
"__=ax.hist(np.log10(iso_cens['halo_mvir_host_halo']), bins=logm_bins, normed=True, \n",
" alpha=0.8, label=r'${\\rm isolated\\ centrals}$', color=mblue)\n",
"\n",
"title = ax.set_title(r'${10 < \\log M_{\\ast} < 10.25}$')\n",
"xlabel = ax.set_xlabel(r'${\\rm M_{\\rm halo}}$')\n",
"legend = ax.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compare 3d clustering of isolated and non-isolated centrals"
]
},
{
"cell_type": "code",
"execution_count": 198,
"metadata": {},
"outputs": [],
"source": [
"from halotools.mock_observables import return_xyz_formatted_array, tpcf\n",
"\n",
"pos_all_cens = return_xyz_formatted_array(\n",
" all_cens['x'], all_cens['y'], all_cens['z'])\n",
"pos_all_sats = return_xyz_formatted_array(\n",
" all_sats['x'], all_sats['y'], all_sats['z'])\n",
"pos_iso_cens = return_xyz_formatted_array(\n",
" iso_cens['x'], iso_cens['y'], iso_cens['z'])\n",
"pos_noniso_cens = return_xyz_formatted_array(\n",
" noniso_cens['x'], noniso_cens['y'], noniso_cens['z'])\n",
"\n",
"rbins = np.logspace(-0.5, 1.5, 25)\n",
"rmids = 0.5*(rbins[:-1] + rbins[1:])\n",
"xi_all_sats = tpcf(pos_all_sats, rbins, period=model.mock.Lbox)\n",
"xi_all_cens = tpcf(pos_all_cens, rbins, period=model.mock.Lbox)\n",
"xi_iso_cens = tpcf(pos_iso_cens, rbins, period=model.mock.Lbox)\n",
"xi_noniso_cens = tpcf(pos_noniso_cens, rbins, period=model.mock.Lbox)\n"
]
},
{
"cell_type": "code",
"execution_count": 202,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEjCAYAAAD6yJxTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsnXtcVNX6/9+bm6IiA+IF8DpYeanMUTqmaR0dMu9pg5fUSguwi546375y/H3P9+g533OOUaey7AYUmqaF4v1SCpaVlqmgXSS7gJqKpDIMIMp9/f7YMzQio1xmGGZY79drXjIze+/17DW+1mc/61nreRQhBBKJRCKR2AMPZxsgkUgkEvdBiopEIpFI7IYUFYnECkVRdIqixCmKstDZtkh+R1EUraIoemfbIbkxUlQkTsM8gK+38V20oigG86vJBnghRAZwCIhoivbMApalKIpQFEVzneO0iqLkm4+Nu96xTYW9fz/zOdGKosSbX9b3qAPWm/spX1GUVEVRdPa4D4l98XK2AZKWh3kwmGZ+q63l+2gAIUSK+b1WUZR4IURME5mY0UTtIISIVRQlD4hD7QtbbesBDTDKLHxOwxG/n6Io0UKIhBrXSAfCLJ8JIQIURdEIIUz2uROJI5CeiqTJEUJkCCFigWQbh8RYDzBCiGzUQbVRmJ+CUxt7HQeQbX5dM0ADKIpiAIxQ7Uk5hLr2j71/v9q8LvP5gTWnvKSgNH+kqEiaFeYBprZpDVND5tQVRdEoirLQPFgahRBNMq1VV8z3lIEqKuG1fK+1+i7NAe3btX8a+PtpgZrTXXAdoZU0X+T0l6S5oQVqexo1og5WdRpYzQPUIvP14oUQLzTGKPN0jNHKxgTLU7NVW4fM32Wbj4uow5SdVgiRpiiKrQFUJ4RIMccu4htzD9bYu3+sqPfvJ4TIUBRlUC1eiHVfUkOUdFj9BpLmgxQVSXMjkN8Hb2tMQIcbnWx+so81X2epPaaLFEWJB+LM0ziWAXk9vwfzE1EH5jTz91lCiDAgpR7NZAGDa7SrB9LM7Wmxg6fiiP6pQYN+v5p2mKf8si19ijnWZPUbZHP1byBpJkhRkbgF5sEyDnVAqxYAO1xXBwy2vp4QwqQoSrZVcFkPRNU870YDtlksLAPwVZ6KZSrI3JbB/HeDBcBR/eMIrLyoUZbPatorhMg2LwC4YT9LmhYZU5E0RwJr+UwD5NXhXBO1Pyk3lMFYTcFYkQUMMv9t5GqbA6l9Cqgmen73PjJQ77H6O6un9AjsF0+xd//URmN+P1DFL7IOU1smanh3EucjRUXS3DjM1YOrhUCus9RXCJEthIhEjTvEmfc5ODrIaxk84zGvbjK3ebiOnkCgZeC0mtaxbPKzFhE90KhVa03YPw36/SyY97Rc5UmZ+6S2JIVGHC+QknoiRUXSrDAPstm1rATSWD25X+/8bHNwPBaIURRlfSN3Yh+m9tVMYfw+0GcDRvM0lb4RK6hMlrZqLAKoNZ5iHmwN9WnAAf1T8/oN/v3MiyFSagiKHlU4alvwMJgm3FMkqRtSVCTOpLZpElCnPxZZ3pjjGvWa/hFCmMx7KaIAnXkHdr0GYPN1MoAM693b5gFzsNVejHAhRIr5lVDrhWpgYzd4NjCtxuCrB0y1BLIXWv1tqK/XYaf+adDvZxbD9dbCYxaPw1Yem8YidrVNg5kFaF1zjg21VBSZ+l7S1JgHwBjUAVMHJADpteyozkadStHaY8mr+ZqR1/MkrALaBiDW0q55ELcMYDWXFOtRVyJZvs8Gki07ymtpIx6YivoEHmu18zwOdRVZttUqramofZBgPtbSpmVXu2VJsD1Wht2wf8zHNer3s+qvQVb3mmWjuQCre16I6s1ZFjHYaxm0xI5IUZFIGoF5QDRgFhnz03cgqiBYvAFHtBuHGvi2LNNdKvdsSJoDUlQkkkZg7V3U8l2qo3bwW0136YF1IFOYSJoHMqYikTSOQ9SS18o8PeWwILJZxIyocQiTFBRJc0F6KhJJIzHHCKxTimhQlwvXKWgvkbgTUlQkEolEYjfk9JdEIpFI7IYUFYlEIpHYDSkqEolEIrEbUlQkEolEYjdafOr7oKAg0bNnT2ebIZFIJC5Fenr6RSFEx5qft3hR6dmzJ4cPH3a2GRKJROJSKIpyqrbP5fSXRCKRSOyGFBWJRCKR2I0WLyo5OTkoisKSJUucbYpEIpG4PC0+phISEkJOTo6zzZBIJBK3oMV7KhKJRCKxH27rqZir2JmACNSCSbLsqEQikTgYt/RULJXpzNXwDmFV2lQikUgkjsMtRcVcayLS/DYcSHWiORKJRNJiaPbTX+ZprPDayrJa1Q0PBLCuX2Eu7WoANLKuhUQikTQNzdZTURRFbxaNGNSiRzW/jwOyhRApZtEIM4tINUKIFCDVfKxEIpFIHEyzFRUhRJoQ4gVsl2SNNouGhVRUAUJRFJ25Gh/m8xc6zlKJRCKRWGi2onI9zPW/a2Lk91rhgwHLMRp+L/MqkbQY0tLSCAsLIyYmptb3zqQutkVGRpKQIGeuXQ2XFBXUGIqxxmcmAEVRLDEUS0wlht+D9hJJi0Gv1xMbG2vzvb2pjwDUxbaYmBj0ev1Vn0mRaf40+0C9DTSYg/NWWEQmEDBZBedTcABVly+jtGqF4unpiMtLJC6FyWQiKyvLrtesKSiOaENif1zVUzHV8plFZGp6MNfFkvvL8qprDrCLiYn8MkrPheWvU37uXH2alEjcCpPJRFRUlN2vmZGRQXZ2tsPakDgGV/VUjFy7IkwD6lLi+lyoobm/2uh0lHz3PRfffJOLb71Fu+HD0UybSrsRI1C8XLVbJa6I9eCbmppKXFwcGs01CybrzAsvvIBOp8NkMpGamkp8fPx120lLS6sWAcu5Fi8jOzub+Ph4wsPDOXToEIsWLaqTbUajkdjYWLRaLfHx8Q1uw9a9SByHS45+QogMRVFqikcgkFbfa1k8lcWLF9crU3G74cNpN3w4ZWfOYEpJoWDDRs48+RRenTrh/+AUNA8a8OkaWl9zJM2A3H//m9Ifjjul7VZ9+9Dl//2/ep0TGRnJ+vXr0el0GI1Gli5dSlxcw1bRJyQkXDNg36gdg8FAdnY2eXl5LFx49ULLiIgI0tPT0Wg0aLVaoqKiWL9+/Q3t0Gq1REZGkp6eDtCgNq53LxLH4arTXwAJNfalRAD1fgwJCQlBCNHg1Pc+XbvS6Zln6P3JHrq+vpxWffuQ93Y8WRER/BoVTWFqKqK8vEHXlkjqQmpqKjqduthx8ODBZGQ0PM2dVqslJiaGhIQETCYT0dHRDW4nISEBrVZb7TXodDrS0ur93NfgNq53LxLH0Ww9FfOyYT1gAAIVRckC0iyJIYUQsYqiLDQLixbIqrFvpU401FO5xl5vb/z0evz0esrPnsW0YSOmDRs4O38Bnh2D0Ex5EE2kAZ+uXRvchqRpqK+n4Gy0Wi0pKSkYjUZMJhNGY73Cileh1+uJi4sjPj6emJgYoqOjq6eM6tuOJahuLSRTp05tsG31beN69yJxHM3WUxFCZAghXhBChAkhAsx/Z9Q45gXzjvoXGpqKpbGeSm14h4bSccF8eu9Jo+tbb+J7623kJSaSpY/g17lzMW3aTGVRkd3ak7RsBg0ahFarJTo6+poVU/UlLS0Ng8FAamoq+fn5ZGdnV08b1bWdlBT12S48PByNRoNer69+2WtQr0sb17sXieNotqLiDiheXvj98Y90e+tNen+yh6D5T1N26lfOLVrEz8Pu5sz8+RR+/DFVJSXONlXiolgC2JZpKevVUg2ZBktNTa0+T6PRVF/3Ru1otVpMpqvDnJY4iPXnjZn+qm8btu5F4lg8W3oZ3SVLlix57rnnALj33nsd1o5nu3a0vfNOAh6eTbvhd6P4tOLSvn0UbNxE/qrVlP6SheLjjU9IiNz7IqkzWq2Wffv2YTKZyM/PJyQkhGPHjpGbm4tWq+XVV18lPT2d0NBQSkpKrnrfr1+/a6534sQJSkpKyMzMJDMzE19fX/R6/XXbMRgM9OvXj1WrVmEymQgODkar1QJqEH3ZsmUUFBSQmZlJSEgIwcHBZGRkXNe2kpISFi9eTGZmJv3790er1da7DVv3IrEPf//7388tWbLkmhkiRQjhDHuaDYMHDxaHDx92StuispLLhw5RuGMHhbtTqSoowNPfH7/77qP9uHG0CR8sBUYikTRLFEVJF0IMvubzli4qISEh4ty5c40O1DcWUVbGpf37Kdz5EUV79iAuX8arY0f8xtyP/7hxtL79dhRFcZp9EolEYo0UFRs401OxRdWVK1zau5fCnTu5tPczRHk5rfv3p0NUFH4Reum9SCQSpyNFxQbNUVSsqSwqonDnRxiTkig7dQqfHj0IfPwx/CdNwsPHx9nmSSSSFootUWnxq78s+1Sa64IFTz8/AqZNRbtzB6HLluHRti25//s3svQR5CWtoPJSsbNNlEgkkmqkp9LMPZWaCCEo/vJL8hLf4fKBA3i0b0/AzIcInD0br8CaiZslEonEMcjpLxu4mqhYc+Xbb8lLfIeitDSUVq3QGAx0mPMo3qEy55hEInEsUlRs4MqiYqE0K4u8d5Mo2LoVhMB//DgCH3uM1jff7GzTJBKJmyJjKjZo7jGVutAqLIyQf/+L3qm7CZw1i8LdqZyYOIkzf3qGslOnnG2eRCJpQUhPxQ08lZpU5OeTv/p98lauRJSXE/jQDDrMm4dXQICzTZNIJG6C9FRaEF4BAXRcMJ+wjz9C88ADGFe/T9bo+8lLWkFVWZmzzZNIJG6MFBU3xrtTJ4L/7x/02rwJ3zsGcP6FF8geO47CnTtp6R6qpOmJjIwkIaFBycQBNVFkWFgYMTExdrSq+bRnb5xlvxSVFkDrm2+me0IC3d59B4+2bTn75//i5PTpXG5EMSeJpL7ExMQ0KqGjXq8nNja2Xuc0RsQa0p69aIzdFpxlf4sXFXcI1NeVdsOG0WvjBoL/9S8qcs5x6qGZnFnwJxnMlzQJlmzHTYXJZKou4uVKuKrdFlq8qDiiSFdzRvH0RPPgFMJ2fUzQ/Ke5tG8fWeMnkPvvf1ORn+9s8yRuiqXuSlMVyTKZTERFRTVJW/bEVe22ptmWE5Y4Fo82bej41FNoIiO5uPx18t9fQ8GmzQTNm0fArJl4tGrlbBOdx0d/gdzvnNN2l9tgzPN1PjwtLY3Y2Fj0ej3h4eGAWpwqNja22iswmUzVtdyNRiNarRa9Xl+nc+2F0Wisvq6lMuMLL7yATqfDZDKRmppa/bkte2vDWqhSU1OJi4tDo9FUFxXLyMiobsdyjezsbOLj4wkPD+fQoUMsWrSousa9yWRi6dKl1f1RF4/B+hyNRkNgYOBVxcxqa8tW399222212p2WlkZMTAwxMTFoNBri4+PZs2cPGo3GZh/Uhq0+tydSVFo4lmB+wOxZnH/xP5x/8UWMa96n44IF+E+YIDMiN3P0ej0xMTHEx8cTFxcHqAO49ftRo0aRnp5efU5kZCSBgYF1OtdeaLVaIiMjq+1ISEi4ZqC3YMve2io3RkZGsn79enQ6HUajkaVLlxIXF1ddETIvL4+FCxdedU5ERATp6eloNBq0Wi1RUVGsX7++um3LYA1w6NChG97bqFGjWL9+PVqtloyMDCIjI6vFyFZbtvo+KyuLiIiIa+y2HJ+cnEx6ejqBVimZbPVBTa7X5/ZEiooEMAfzExMoPnCA8y/+h3N/WYRxxUo6PfdftL377pZVy6UenkJzwdqzCAwMrB6UU1JSrvE6pk2bxtKlS6sHUlvnOtremJgYYmNjmTp1KtHR0XW215rU1NTq4wcPHlzrMdZYPCCLaOh0uuryw5Z/rZ/yw8LCrtsflnMsNuh0OlJTU2/YlnU/WLD0fYcOHWy2ZzneYDBUf1bXPrDV5/ZGiorkKtoOGULP9eso/OgjLryyjNNR0bQZMoROzz2H7639nW2exAaBNZKJGo1GQH3SrvmdZcrkRufWRkJCwg1Fpy7TZ3q9nri4OOLj44mJiSE6Opr4+Pg62WuNVqslJSUFo9GIyWS6ru3w+3SW9eA+depUQJ1Kq9n2jajtHMu9X68tC7X1fVhYmM32auvXuvaBrT63N1JUJNegeHjgP24c7SMiyP8wmYtvvcVJg4H2Y8fS8Zk/4dO9u7NNlNSRsLCw6idnCyaTqdappLpgr6fbtLQ0DAYDBoMBk8lEZGQk2dnZ9bZ30KBBJCYmYjAYyMjIIDk5udbjUlJSMBgMhIeHk52dfVWMxvK3Tqezeb4trnfO9dqqKxa7r0dd+8BWn9s7ftbiV3+1pCXF9UXx8SHw4dmEpe6mw7wYij75hKxx48n957+ouMEToaR5EB0dfc3ceXJyMosWLXKSRSqpqanV3odGo6kWjfrYawnGWwfF4feVZlqtFpPJdNU5lliL9ecWT8Iy4Ft/dyOvTK/XYzKZrrLZssfkem1dj9rstsWN+sAaW31ubzxb+mD6wQcfLMnJyeHee+91tinNFg8fH9oOGYL/5ClUFRVhWrcO0wcfQFUlrfv3R/H2draJLZaMjAwWL15MZmYmoaGhlJSU8Pzzz7Nv3z569+5Nv379iIiI4Pnnn6egoIDt27czYcIEhg0bVqdzHWFn//79URSFkpISMjMzyczMxNfXt3pQv569r776Kunp6YSGhjJ+/Hj27duHyWQiPz+fkJAQjh07Rm5uLgaDgX79+rFq1SpMJhPBwcHVT+QREREsW7aMgoICMjMzCQkJITg4+KrvFEWpXlW1bt06AgMDGTRoUK339sgjj7BkyRJKS0vJzMxk2LBhBJjz7Nlq63p9P3HiRA4fPnyV3WlpaSxbtoyMjAx8fX2rbdFqtTb7QKvVXtVfvr6+Nvu8Ifz9738/t2TJkmt2acqEkm6YUNLRlGZnc/7ll7mUtgevjh0JWjAfzYMPoni0eMdXImkxtLiEkoqiRJtf8YqiNN023hZAK62Wbq+/To+1a/Du2pXc//0bp594Qk6JSSQS9xQVRVF0wGEhRAKw3vyS2Jk2Oh091q6h89/+l8tfHeDEpAcoPvC1s82SSCROxC1FBdACltSch83vJQ5AURQCH3qInuuS8WjXjl/nzOH8q68iKiqcbZpEInECzV5UFEUxKIpS6/ZeRVEWmr+PVhSleq2jECIFsKTn1AM3XnIhaRSt+/Sh14YU/KdMJu+ttzn18COU5+Q42yyJRNLENFtRURRFryjKQlSP45pENmahyRZCpJinucIURale0C2EsKzJmwa4doY2F8GjTRtC/vUvQl58kdIffyR78hQKa+w5kEgk7k2zFRUhRJoQ4gXAVtGPaLNHYiGV36e8ANWTAaKsBEbSBPhPGE+vTRvx6daNs/MXkPuPf1BVWupssyQSSRPQbEXlepgD8TUxok51WY4xAAlCCJOiKA1fjC1pED7du9Nz7RoC58whf+0HnJw6jVIXrhEhkUjqhkuKChCIKiLWmAAURdGYRScRSFcUJZ/f4yuSJkTx8aFz7EK6JcRTcf48JwyRmDZskKWMm4ia5WRdvTwuuMc9uDuuKioaVGGxxiIygUKIDCFEgBAizPxvRBPbJ7Gi3YgR9Nq8Gd8BAzj3P38l57+eo7KoyNlmuT01y8k6szwuuHaJXEndcVVRqS1GYhEZuQOvGeLduRPd332Hjs88Q+GuXZyYPIWSzExnmyVpIly9RK6k7riqqBi5dkWYBq5a9VUnLAklLa+WngvNkSiengTNi6HH6tWIigpOPjSTgu07nG2WxMG4Q4lcSd1xydT3QogMRVFqikcgjdiPsnjxYikoTUQb3UB6paznzJ+eIee55yg9/gMdn31WVplsIPUpJ1sX7FkeNzY2loyMDJcvkSupOy4pKmYSFEUxWC0rjgDq/b8pJCSEHLlJr8nxCgqix4okcpcuJe+ddyn54TihL7+Ep7+/s01zOepaTrau2LM8ruV9baV9XalErqTuNFtRMa/g0gMGIFBRlCwgTQiRASCEiLXsqEdNw5JVY9+KpJmj+PgQvHgxrfv2Jff//smJyKl0e+N1Wt10k1PteuaZZzh69KhT2r7jjjtYtmxZvc6pb0nd6+GI8rg3whVK5ErqTrMVFbN4ZAAvXOcYm9/VFUtMRU5/OY+AqVNp1bs3Zxb8iZPTphPyQhx+jajz0NKob0nd6+GI8rg3whVK5ErqTrMVlaZCTn81D9rodPTakMKZ+Qs48/R8gp56iqCnnnRKjZb6egrOpq7lZOuCo8vjgmuWyJXUHVdd/WU3ZDnh5oN35870WL0K/8mTufjGG5x5ej6Vly4526xmTX3KydYFR5THBdcvkSupO7KcsCwn3KxQvLxoN2oknv7+5K9ZQ1FaGm2H3oWXuTyr5GrqU062pKTkqve2ygXbuzxu7969MRgM15T2bY4lciV1R5YTtkGDywmnvwcnPoNbDdB7FHi1sr9xLZziA19z9plnEJWVhL70H9qNGOFskyQSiZkWV064rjR4+qu0CLI+hQ9nwIs3weanIOsTqJTFqexF2yF/oGdKCt6hoZyOmcfFhESZN0wiaeZIT6WhngpAZTlk74XvN8AP26GsCNp2hH4PwK0PQrc/gBMCze5G1eXLnPvrXync+RF+Y+4n5J//xKNtW2ebJZG0aGx5KlJUGiMq1pRfgZ9TVYH56WOoKIH2XeHWyarABN8BitL4dlooQgiM777L+ZdfoVVYGF1fX45Pjx7ONksiabFIUbFBSEiIOHfunH33qZQWwY8fwXcpkLUHqiogMEwVl36ToFNf8JApSRrCpf37yfnzfyGEIPTFF2h3zz3ONkkiaZFIUbGB3TwVW1w2wg9bVQ/mxBeAAM9WEHQzdLwFOvaBTn3UfwN6gWeL3zp0Q8rOnOHM/AWUHj9O0PynCZo3zyn7WSSSlowUFRs4XFSsKcqFX/bAhR/g/HG48CMU/Pr7954+0KG3KjAd+6ii06kvBGrB07tpbHQRqq5c4dzfFlO4bRvtRo0iJO55PNu1c7ZZEkmLQYqKDZpUVGqj9BJc/EkVmAvHf3/lnwLMv42HNwQPgO5DoPtd6r9tg5xnczNBCEH+6tX8FvcCPt270/WN12kld1JLJE2CFBUbOCSmYg/KLkPez6rY/HYMTh+Es+lQWap+3+Gmq0UmUNtiFwIUf32Qs88+iygtJSTueZk3TCJpAqSo2MDpnkp9qCiFnKPw61fw6wE4fQCu5Kvfte10tch0ub1FxWfKz53jzII/UfLdd3SYF0PH+fNlfRaJxIHIzY/ugFcr6P4HuPsZeOhD+O9sePJrGL8MwkbCuW9g1yJI/CM83x3WPwon90MLeHDwDg6mx/ur8X9wCnlvx3P6iSeoLChwtllNQmRkZKPqv6elpREWFkZMTIwdrWo+7dkbV7ff0UhRcWU8PNSVY4PnwJR4eOZb+PMPYFgBA6arO/5XjoW3hsHhJDV+48Z4tGpF8D//SZcliyn+6gAnIqdS8uNPzjbL4cTExDQq35Veryc2NrZe5zRGxBrSnr1ojN0WnGm/KyBFxd1oHwK3ToHxL6sCM/F1dU/M9mfh5b7w0V/g4i/OttJhKIpCwPTp9HjvPcSVK5ycPp3Cjz5ytlkORa/XN2mqd5PJVF1bxZVwVbtdjRYvKm6d+t6nDehmQ8zn8Fgq3DwaDr0Drw+C1ZPh+E6oqnS2lQ6hjW4gPTek0LpPH84++2fy3k1yy7xhlnTwTVVG12QyERUV1SRt2RNXtdsVaTmRXBu0iCJdigLd7lRfo/+tZlg+nKQmw/TvDuFzYeDD0LaDsy21K96dOtF95QpyYv/C+RdfpDw3l85/iXWrAL7RaCQ2NhatVltd8fCFF15Ap9NhMplITU2t/txkMlWXBDYajWi1WpvTZtZClZqaSlxcHBqNprrWSUZGRnU71vXh4+PjCQ8P59ChQyxatKi69LDJZGLp0qWEh4cD1MljsD5Ho9EQGBh4VY2V2tpKS0sjNjYWvV5f3VZqaiq33XZbrXanpaURExNDTEwMGo2G+Ph49uzZg0ajsdkHtWGrz1skQogW/Ro0aJBokVSUC3FssxArxgmxuL0Q/+goxMZ5QmR/LkRpsbOtsytVlZUi99//Fpm39BGn5y8QlSUlzjbJrsTHx4vo6Ojqv1NTU6u/i4uLq/5bp9NddZ7BYBDp6enXXEMIIbRa7VXfLVy48KprWr+3Pic/P18IIUR6erowGAxXtW35TgghFi5ceFV7taHT6URWVlb19bRabZ3aio+Pv+peLfbbsjsuLq76+PXr11df93p9UNc+d2eAw6KWMbXFeyotFk8vNQ9Zv0lw/gd1WuybD+GbtaB4Quf+0DXc/Bqs5i5z0VQoiocHnRctwis4mPPPx/FrXh7d3ngdTxtPnX/fdozMnMImtlKlX0h7Fk/o3+DztVotMTExxMbGMnXqVKKjowG1hG/NuMu0adNYunQp69evv+Y6qamp1ccPHjy41mOssXhAlid5nU5XXRXS8q/1U35YWBjp6ek2r2c5x2KDTqcjNTX1hm1Z94OFwMBA0tPT6dDBtiduOd66zHFd+8BWn7dUpKg0kPz8fFq3bo2vr6+zTWk8nfrCuJdg1GI4tR/OHIYzh+DbdXD4XfWY1hoIHfS70ITqoE2gc+2uJx0efRTvzp3JWRjLyZmz6J4Qj3doqLPNsit6vZ64uDji4+OJiYkhOjqa+Ph4Dh06RGDg1b+XZYqnNrRaLSkpKRiNRkwmE0aj8brtWqazrAf3qVOnAupUWs22b0Rt51gG+Ou1ZaHmuUajkbCwMJvt1bbQoa59YKvPWypSVBrI0qVLeeutt5g0aRJTp05l9OjRtGrl4tUfW7eHW8aoL1CD+Bd/UgXmzGH19Vkc1eljOvSG0MHQLRx6joCgm5r9rv72Y8bg2aEDZ556mpPTZ9AtIZ7WfftedUxjPAVnk5aWhsFgwGAwYDKZiIyMJDs7m7CwsOonfQvWdeBrMmjQIBITEzEYDGRkZJCcnFzrcSkpKRgMBsLDw8nOzr4qRmP5W6fT2TzfFtc753pt1RWL3dejrn1gq8+bckVec8I15zOaARMnTmT69Ol89NFHTJo0iU6dOvHII4+wc+dOysrKnG2effDwVL0Y3cMw8TV48ktYdBqL7kFKAAAgAElEQVQe2Qaj/gZBt6jVLnf8F7wRDi/1gQ1RcOR9MP164+s7ibZ33knPtWvA05NTs2Zzaf9+Z5tkN1JTU6u9D41GUy0a0dHR16wQS05OZtGiRddcwxKMtw6Kw+8rzbRaLSaT6apzDAYD2dnZV31u8SQsA771d9eb+rKcYzKZrrLZssfkem1dj9rstsWN+sAaW33eUvF0y6W09WDJkiVLnnvuOQDuvffeOp/XvXt3Jk6cyJ///GeGDh1KVVUVmzdvJikpiTfeeIOffvqJ1q1b06NHDzzdaLURXq0goAf0GAq3PQhD58PtU6FzP0DAic/gu/Vw4C01RnM+U81j1jYIfJpPtUavDh1oP2YMl/Z+hnH1anxCQ2ndp4+zzao3GRkZLF68mMzMTPr374+iKJSUlJCZmUlmZia+vr7Vg3pERATPP/88BQUFbN++nQkTJjBs2DAyMjJ49dVXSU9PJzQ0lPHjx7Nv3z5MJhP5+fmEhIRw7NgxcnNzMRgM9OvXj1WrVmEymQgODq5+Io+IiGDZsmUUFBSQmZlJSEgIwcHBV32nKEr1qqp169YRGBjIoEGDar23Rx55hCVLllBaWkpmZibDhg0jICDgum1Z90doaCglJSU8//zz7Nu3j4kTJ3L48OGr7E5LS2PZsmVkZGTg6+tbbYtWq7XZB1qt9qr+8vX1tdnn7szf//73c0uWLLlmN6nM/WXH3F+lpaWkpqaybt06Nm/eTFFRER06dGDKlClMmzaNe+65By8vN59xFEIN/J/4XH2d3Ael5nQpHftCrxHqq+fd4Ft7oLwpqSwq4szT87n89dd0fPZZOkRHoTTzKTyJpDkgE0rawFEJJUtKSti1axfJycls3bqV4uJiOnXqxOTJkxk5ciTDhw+vfopza6oq4dzR30Xm1FdQcQUUD+g5HG4zQN8J4BvgPBPLyji36P9RuGMHmhnT6fLXv7rVXhaJxBG0SFFRFEULGIUQNidSmyJL8ZUrV/joo49ITk5mx44dFBcXA+qyyuHDhzNixAiGDx9OWFiY+z8lV5SqKfx/SYPvN0L+CbVezE0RarnlW8Y4ZZpMVFVx4eWXyXvnXdqNGkXof17Ewx1W9kkkDqLFiYqiKHogFogTQtiM4jV16vvy8nKOHDnCF198wRdffMG+ffvIy8sDoEuXLgwfPrz6ddttt7lXPKYmQkDOEbXU8vcboSgHvNuownKrAXqPUmM4TYjx/TX89q9/4Xv77YQufw3vTp2atH2JxFVocaICoChKPLC+OYlKTaqqqjh+/Diff/55tdCcPn0aAH9/f4YNG8bw4cO5//77GTBggPt6MlVV8OuXqsAc2wxXjNDaH/pOVKfIeg5XV6M1AYWpqeQsjMWzXTtCX32VNrqBTdKuROJKuKyoKIpiAMKFENfkmlYUZSGQDQQCCCESanzf7EWlNk6dOlUtMJ9//jnHjx8Hfl9xNnHiRO655x58fHycbKmDqCyH7L3wXQoc3w5ll9QiZP0nq2n+O/W94SUaS8mPP3Fm/nzKz52j86K/EDBjhvsKukTSAFxOVMzTVzogAsgWQsTU+D4OOCSESKntvfkzlxSVmuTm5rJz5062bt3K7t27uXLlCn5+fowZM4aJEycyZsyYeu9YdhnKr8BPu1QP5qddUFmmpva/5y/Q8WaHNl1ZWEjOfy/k0mef4T95Ml0W/w2P1q0d2qZE4iq4nKhYMIuFphZRyRdCBFi91wOxQogIq8/cQlSsuXLlCnv27GHr1q1s27aN3NxcPD09GT58OBMnTmTChAn07t3b2WY6hstG+HI5fB2vriC7bSrcsxA62E6/0VhEVRUX33iTi2+8Qet+/ei6/DW3S+0ikTQEtxIVRVF0wJ4aoqID0oUQitVnbicq1lRVVXH48GG2bt3K1q1b+e677wDo27cvEydOJDIyEp1O537TNsUXYf+rcDBR9VwGzIARz0FgL4c1WfTJp+QsXIji5UXoyy/RduhQh7UlkbgC7lajPhComd3NBKAoisb8rwEYDESaBcft8PDw4M477+Sf//wn3377LSdOnOC1114jNDSUl156icGDBzNgwABeeeUVLly44Gxz7UfbILjv/+BP38Af5sH3KfD6YNi6wGHpYfxG/pFeKevx6hjEr49HkffOO25Z9EsiaSz19lQURRmJGufQYA6Qow7wWUCaEOKoXQ2s3VMxAIk1PBUNkA+ECSHqXAbPlT2V65Gfn09ycjIrVqzg4MGDeHl5MX78eObMmcOYMWPw9vZ2ton2o/Ac7HsZ0leqy5QHPQJ3/xn87T9NVVVcTM7//JWijz/Gb/Rogv/1LzzbNZ/0MxJJU9Go6S9FUfyBOKAXkAqcQF11ZdlUqAG05lc4kIca32h0UQoboqJHndayFhUtqrAFXG+zY01CQkLEuXPnqt8vXrzY7UoLHzt2jJUrV7J69Wp+++03OnXqxOzZs5kzZw79+7tuRt5rKDgDX7wEGavVHfuD58Ddz4JfF7s2I4TAmJTE+Zdexkfbi67Ll9Oql+Om3iSS5kiDRUVRlAdRxSRRCFFQx8b8galAnhBiYwPstb6WrZhKzfjJNZ/VBXf1VGqjvLycjz/+mBUrVrBt2zYqKioIDw9nzpw5TJ8+vTpZn8uTfwq++A8cWQOe3jBoDgx5Qk2EaUeKv/qKs8/+GVFRQcgLcfiNHGnX60skzZkGxVTMgpIhhPhPXQUFQAhRIIRIBPYoijKl/ube8PoZ/O4lWQgEbpz/ugY5OTkoiuJ23klteHt7M2HCBDZu3EhOTg6vvPIKpaWlPPnkkwQHBzN9+nR27txJaWmps01tHAE9YOJymH9YTf1yKBFeuwPWPaLWhLETbe+6i14bUvDp0YMzTz7FhddeQ1RU2O36Eokr4pKrv6w+v+4+lbrQkjyV2hBCcOTIEVasWMHatWsxGo34+fkxduxYHnjgAcaMGYO/v7+zzWwcBWfhYDwcXqlmTO42BO56CvqMs8su/aqSEnL//g8KNm3Cp1cvOv75Wfz0evdbdSeRWOFyS4rN01l6IAbVC1mKuhAgw+oYy456LWCquaO+LlhiKu4YS6kvpaWl7Nmzh82bN7NlyxbOnz+Pt7c3I0eO5IEHHmDSpEmunVm59JJaQOzAm2A6BQG9YMiTMHBmo5NYCiG49MknnH/pZcqys/EdOJBO//0cbVp4wSaJ+2I3UTGv/tIKId6xl3HOpKV7KraorKzkwIEDbN68mU2bNlXXBR8yZAgPPPAADzzwALfccouTrWwgVZXwwzb46nW1VHJrDQyeC3dGQ/vGiaaoqMC0aRMXl79OxfnztBs5kk5/fpZW7rohVdJisaeoPA8MFEKMtpdxzkR6KjdGCEFmZma1wFhKwfbp06fagwkPD3fNjMqnD6q79I9vB8VTTV5511PQ5bZGXbbqyhWM760i7513qLp8Gc2DUwh6+mm8O3e2k+ESiXOxp6g8KITYYDfLnIz0VOrP6dOn2bJlC5s3b2bv3r1UVlYSEBDAyJEj0ev1REREEBbmuNQpDsF4Qi2BfOR9KC+GXveoK8ZuGg0eDd8jXJGfT97bb2Nc+wGKpyeBDz9Mh6jH8fTzs6PxEknTY09RGQj0auxS4eaCFJXGYTQa2bVrF6mpqaSmpnLmzBkAevXqRUREBBEREYwcOdJ1El5eyVc3UX6doNZ3CegFf4iBO2ZC6/YNvmzZmTNcWPYqhdu346nREPTEPDQzZuDhrpmmJW6PXT0VIBF1o2EyDthF35TI6S/7IYTgp59+qhaYTz/9lKKiIhRFYdCgQdUiM3ToUFq1atriW/WmslyNu3z9Npz+GnzaqcJyZzQENTw+UpKZyfn/vETxl1/iHRpKx2f+RPtx41Aa4Q1JJM7AnqLy36j7QQajpmvRA/7mz9a7WgBfeiqOo7y8nEOHDlWLzIEDB6isrKRNmzbcc889TJ8+HYPBQJs2bZxt6vU5m6FmRv5+A1SVw033qd5L2Cho4LLhS/v3c/6llyjN/IFWffrQ8Zk/0e6ee+QyZInLYE9RiQKyhBCfWH2mBUah5t36S2ONbUqkqDQdhYWF7N27l7S0NHbs2EF2djZ+fn7MmDGDuXPncueddzbvQbXoN0hfAYfeheLzEHSz6rkMmAGt2tX7cqKqisIdO7mwfDnlv/6Kr05Hp2efoU14uAOMl0jsi133qSiK0gs1rvLJDQ9u5sjpL+cghOCLL74gKSmJ9evXc/nyZfr378/cuXOZNWsWnZpzbfiKMji2Cb5+C3KOQCt/0M2GO6MgoGe9LyfKyzFt3MTFN96g4vx52t59Nx2feQbfW90oL5vE7WiQqCiK0lMIcbKRDTf6Go5EeirOp7CwkOTkZJKSkjhw4ABeXl5MmDCBxx57jNGjR+Pl5eVsE2tHCHWfy4G3IHMLINS0MMOfg0596n25qpIS8td+QF5CApUmE36jR9NxwXxaudpKOkmLoKGi0qiVXuagflZzDuRLUWleZGZmsmLFClatWsX58+cJDg7mkUceYc6cOdx8s2PLBzeKgrNqUP/Qu1B+Gfo/ACP+GzrX39uovHQJ44qVGFesoKqkBP9Jk+j49FOy4qSkWdGYLMX+wCLgF2DdjdLZK4rSHpiGmjolvjl7KSBFpblSXl7Ojh07SEpKYufOnVRWVjJ8+HCmTZvG2LFj6dVcU80X58GBN9QlyWVF0HcCjFgIwbfX+1IVRiN5CYnkr12LEIKAadMImheDV1CQAwyXSOpHo2MqiqKMAiJR83AJ1JxbeeavO6CKSABqoax4V4m3yJhK8ycnJ4fVq1ezYsUKfvzxR0DdzT927FjGjh3L8OHD8Wlu+z0uG1XP5cDbahLLW8aqnkto/XOBlZ87x8U338K0cSOKj4+6gXLuHDxdPdGnxKVxRKBei1qcC9Q09NlCiBONstIJSE/Ftfj555/ZuXMnO3fuZO/evZSVldGuXTv0ej1jx45lzJgxdO3a1dlm/s4VExxMgK/egBKTuhx5xELoVv8VXqUnTnBx+esU7tyJ0qYNmgcmETBrNq20zdRrk7g1LpeluKmQouK6FBcX88knn7Bz50527NjB6dOnAbj99turvZi77rqreQT6SwrVui5fvg5XjBA2Eu6Jhe5D6n+pH3/EuGIlhTt2IMrLaTt8OIEPz6btsGFyE6WkyWhsOeFmvYKrMUhRcQ8sSS8tXsy+ffuoqKhAo9EwZswYDAYD999/v/M3WpZegsPvwv7X4PJF6DUC7poPvUfVu7ZLxcWL5Ccnk//hh1ReuIhPr14EzJ6FZtIkPNo2LpW/RHIjGisqU2quAFMUpb2toL2iKM+Z/0ywR516RyJFxT0pKCio3mS5bds2Ll68SJs2bRg3bhwGg4GxY8fSrl39NyzajbLL6kbK/a/BpVzw7wa6h2HgLGgfUq9LibIyCnftwrhqNSXffYeHnx+aBx8kYNZMfJrTVKDErWjM6q+BQEDNwLuiKCNrC8abd9y/DcwDdEKIJxpluYORgXr3p6Kigs8//5yUlBQ2btzIb7/9RuvWras9mPHjx9O+fcOTRTbOuDL46SM4vAKyPwXFA26+HwY9Cr319fJehBBcOXqU/NWrKdy1G4Sg3cg/Ejj7YdrcGd68sxVIXI7GiEov1PopNT0VW6JyGLWs7xOukCZfeioti8rKSvbt20dKSgobNmzg3Llz+Pj4MHr0aAwGAxMnTkSj0dz4Qo7AeAIy3lPT7xdfgPZdVe9FN7ve3kt5bi75H3yIKTmZSpOJVrfcQsDMh2g/erRcNSaxC42d/vIH9gAfmP/VAuFCiEU1jhsIHEbNAXZSUZQoIUSiPW7AUUhRablUVVXx1VdfkZKSQkpKCmfOnMHb2xu9Xk9kZCRTpkzB3xkDsMV7SV8JWZ80ynupKimhcPt2jKtWU/rTTyje3rQdMQL/8eNod++9ePj6Ouw2JO6NPfap6FAzEfsDGcDzgLD2YBRFWQf4u1JVSCkqElAF5tChQ9UCc/LkSVq3bs3EiROZPXs2o0ePxtvbu+kNM56AjFVm7+X8797LwFngX/cd9kIISr4/RuH27RTu3EnFhQt4tGlDO/0o/MePp+1dd6E44/4kLos9sxT7CyEKzH8/iJr6PhW4E1iIGkdptmlZaiJFRVITIQRff/0177//Ph9++CF5eXkEBQUxffp0Zs+eTXi4E+ITleXw40dqcN/ivdw0+nfvxbPuy6ZFZSWXDx2mcMd2CnftpqqwEE+NBr/7R+M/fjy+Op1cmiy5IQ7bp2KeGotG3VWfLIQ40qgLNjFSVCTXo6ysjI8//pj333+frVu3Ulpays0338ysWbOYNWuWc9LF5J+EjNVwZDVc+g3ah8LA2ar3oulWr0tVlZVRvG8fhdt3UPTJJ4iSEryCg2k/dgz+48bRqm9fGeCX1Irc/GgDKSqSumIymdiwYQOrV6/ms88+A+Duu+9m1qxZTJ06lYCAgKY1qLIcftqlei+/7FELhvWOUL2Xm+6rl/cCUFVcTNEnn1K4fTuX9u+Higp8wsLQRBrQTJ4sA/ySq5CiYgO5pFjSEH799VfWrFnD6tWr+eGHH/Dx8WH8+PE8/vjj3HfffXh61m8jY6PJP6V6Lkfeh6Jz4Besei+62aDpXu/LVeTnU7RrNwWbN3Pl6FGU1q1pP3YsATNm4HvbrQ64AYmrIUXFBtJTkTQGIQRHjhzh/fffZ82aNZw/f54ePXoQFRXF3LlzCQ4OblqDKivg593qyrGfd6uf9dbDoEfUFWSe9Q/Gl/zwA/kffEjB9u2Iy5dpfdttBMyYQfuxY/Bo3dq+9ktcBikqNpCi4t7s+PYcu47lcnPndtzc2Y9buvjRLaANHh72jxOUlZWxZcsW4uPj2bNnD15eXkycOJGYmBj0ej0eTR38Np1WPZeMVVCUo3ovd0bD4DngW/+pusqiIgq2bCX/gw8oy8rCw98fzeTJBEyfhk/Pnva3X9KsaXGioiiKATV7shZIE0Jk13acFBX3Zu7KQ+z98TxVVv/Nfb09uckiMp39uLmL+m/n9q3sFpT++eefSUxMZMWKFVy8eBGtVktUVBRz5syhc+fOdmmjzlRWwC+p8HW8umvfu606LTbkiYaVPxaCy4cOkf/BBxSlpkFFBW2HDiXgoRm0u/delOaQwFPicFqUqCiKogVihBCx5vfrhRCRtR0rRcW9mbvyEBeKSvkwegg/n7/ET7lF/PhbET/9VsSPuUWcLyqtPrZ9ay9u6eLHzZ390HUP4IGBoXg20qMpLS1l48aNxMfH89lnn+Ht7c0DDzxATEwMf/zjH5vee8n9Tk3D/10KiEq1iNjQBdD1mrGhTlRcuIApJYX85HVU5Obi1aULmqmRtL9/DD69esqVY25MSxOVhYBJCJFgfp8lhKi10LcUFffGIirb5t9d6/f5xWWqwJhFxiI2hSUVDOyu4UXDAHp3sk/iyePHj5OQkMDKlSvJz8+nd+/eREdH8/DDDze991KYo9Z5OZwEJQXQbQgMfVotJlbPbMkAoqKCS599Rv7aDyjevx8A75AQ2t59N22H303bIUPw9POz911InIjLiop5Givc4nXU+G4hagXKQAArEYlDzT+WYn6fBQwSQphqXkOKintzI1GpDSEEW7/JYfHWY1wuq+S/Im7m8eHaRnstFq5cuUJKSgrx8fHs378fDw8P9Ho9M2fOZPLkyfg15eBbekmNuxx4E0ynIFALQ56EO2aCT8PKBJSdOUvxF59zad9+Lh84QFVxMXh64nvHHbQbfjdth91N6/795AZLF8flREVRFD2gAyJQq0rG1Pi+pnBUvzf/nWXtqSBFpUXSEFGxcL6ohL9u+p7dmb9xRzcN/4m8nd6d7DvgZ2ZmsmbNGtauXcvJkyfx9fVl4sSJzJw5k9GjRzddmeTKCji+TS0idvawGsgf/Jga2PdruBclysu5cvQol/btp/iLLyjJzATAMyCAtsOG0fbuYbQbNgyvjh3tdSeSJsLlRMWCWSA0tYhKvhAiwOq9HogVQkTUMv111bHWSFFxbxojKqB6Ldu+PcfiLd9TXFbJs/qbiRreCy9P+z5lCyH48ssvWbNmDevWrSMvL4/AwECmTp3KzJkzGTp0aNPEX4SA01/Dl8vh+A51CXKf8eqS5J4joJE2VOTlUfzllxTv28el/V9SefEiAK369MF/wngCZszAw9mF1CR1wq1ExZzcck8NUdEB6UIIxRyojxVCxCiKogESZaC+ZdJYUbFwoaiUv235no++z2VAV39ejBzAzZ0dM01VVlbG7t27WbNmDVu2bOHKlSv06NGDhx56iJkzZ9K/f3+HtHsNeVlwMBG++QBKTOpKMd3D6tSYX5dGX15UVVF6/DiX9u3n0mefcSU9Hc+gIIJiYtBMm4pHU3lpkgbhbqKiB+Ktg+9mIclCLShmUhQlGjXeogNS5JLilom9RMXCjm/P8b9bvudSSQV/0t9EzAit3b0Wa4qKiti8eTNr1qwhNTWVqqoqBgwYUJ17rEuXxg/uN6S8BH7YptZ6OfkFKJ5wyxjQPdKgMsi2uJyezoVlr3L50CG8goMJemIemsmTZfbkZoq7iYoB1fuw9lQ0QD5qLZdaBaQ2pKi4N/YWFYC8S6X8besxdnx7jttC/Xkx8nb6dHF85cjffvuN5ORk1qxZw8GDB/H09OT+++/nkUceYcKECbRuit3teVmquBxday4kFqomshw4u97JLGtDCMHlAwc4v2wZJd98i3f37nR8+inajxuH0tSpbyTXxZaouOryi2sC7phXgAHG+lwoJycHRVGqXzL/l+RGdGjXijce0vHmTB05pitMWL6P5Xt+pryyyqHtdu7cmQULFvD111/zww8/sHDhQo4ePcrUqVMJCQnhqaee4uDBgzj0QbFDGET8A57NhKmroGMf+OwFWHYbvP8gZG5VE102EEVRaHvXXfT88EO6vv0WHm3bkrMwluyJkyj8eBeiyrF9LGk8ruqpVMdPrvdZXZCeinvjCE/FGmNxGYu3HmPbNzn0D2lP7P19GH5TUJNt+qusrGTPnj289957bNy4kZKSEvr27cujjz7KrFmzCAmpXxniBpF/Sl2WfOR9NR1M206q9zJ4ToOSWVojqqoo2p3KheXLKcvKolW/vnRcsIB299wjN1Y6GbfyVIQQGVzrrQSiVqasFxZPRXoo7omjH5oC2/qwfMZA3p41CGNxGQ8nHWTi6/v5+PtzVFU5/oHN09OT++67jzVr1pCbm0tiYiKBgYHExsbSrVs3xowZQ3JyMiUlJY4zIqAHjPwfeOY7mJEMoYNg/zJYdjusnQY/7YaqygZdWvHwoP39o9Fu3ULIC3FUFV3izLwnODV9BsVffWXnG5HYA5f0VKw+r3WfSn2uLz0V92bOioNcvFTmME/FmtKKSjYfOctbe7M4mXeZ3p3a8cQ9YUy8IwRvBwbza+Pnn39m1apVvPfee5w+fRp/f39mzJjB448/zqBBgxxvgOm0mik5Y5VaBlnTHQbNUWMv7Rq+J0WUl2PatImLb75FRW4uvjodAdOm4nfffXj4+trPfskNcblAvXk6Sw/EoHohS1ETQ2ZYHWPZUa/Fal9KfZD1VNybOSsOkldcxtanHS8qFiqrBDu/O8cbn/7C8dwiQjW+xNyjZergbrT2btpgc1VVFXv37mXlypWkpKRw5coVBg4cyOOPP87MmTPxd3ThrYoyOL5dTQdz8gvw8IZ+kyD8Meh+l1pYrAFUlZZiWrce4/urKT/1Kx7t2tF+/Dg0DxpofWt/OTXWBLicqDQV0lNxb5whKhaEEHz643ne+DSL9FP5BLXz4bG7tcwa0h2/1k2/TNZkMvHBBx+QmJjIkSNH8PX1ZerUqURFRTF06FDHD8QXflTF5egHUFoAnfrB4Llw+zRo3bDVc5aMyQUbNlC4azeipIRWt9yCxmDAf8J4PDUaO9+ExIIUFRtIT8W9caaoWBBCcPCEkTf2ZvH5Txfwa+3Fo0N7MmdYLwLbOmeDX3p6OomJiaxdu5aioiL69u3L448/zsMPP0xQUJBjGy8rhu83wKF34Nw3air+26eqKWE692vwZSuLiijcsQPT+hRKjh1D8fHBT69HY3iQNkOGyFxjdkaKig2kp+LeNAdRsebbMybe/DSLj4/l4uvtyfQ7u/Ho0J706NDWKfZcunSJdevWkZiYyIEDB/Dx8WHy5MlERUU5PjW/EHA2Aw6/q4pMRQncNBrufhZ63NWoS5ccP44pZQMF27ZRVVCAd2go/g9OQTN5Mt5NXY3TTZGiYgMpKu5NcxMVC7+cL+KtvdlsPnqWyirB4B4BTNF1Zdztwfj7OmcH+ffff88777zDqlWryM/PR6vV8thjjzFnzhzHl0W+bFQ9l6/fhst50O0PMOwZtQRyI4StqrSUotQ0TBtSuPzVAVAU2o4YTlDMPNroBtrxBloeUlRsIKe/3JvmKioWcgtK2HTkLBsyzvDL+Uv4eHkQ0a8zD+pCGX5TxyZfNQZQUlLCxo0bSUxMZO/evXh6ejJ+/Hiio6MZPXo0no7c2V52Wd3v8tVyMP2qbq4cugBuiwSvxk0Vlp05Q8HGjeQnr6MyL4+2I4bTccGf8L21iXKpuRlSVGwgPRX3prmLigUhBN+dLWBjxlm2HD1L/uVygtr5MOmOUKboQukX3N4pK5p+/vln3nnnHVauXMn58+fp1q0bc+fOZe7cuXTv3riNjdelsgKObYJ9r8D5Y2o6mLueUvONtWpc0bSqy5cxrlmD8Z13qSwowC9CT9DT82l9y812Mr5lIEXFBlJU3BtXERVryiqq+OynC2xIP8Oe479RXino08WPKbpQHrgjlE7tmyDHV02bysrYtm0biYmJ7N69G4D777+fqKgoxo8fj7ejkj4KAb+kwb5lcGoftNbAnVHwh3nQtnELCiovXcL43nsYV6ykqriY9mPHEvT0U7Tq1ctOxrs3UlRsIKe/3BtXFBVr8ovL2EoJU1cAAB/eSURBVP7dOTakn+HoaRMeCgy/qSPTw7txX/8udqtGWR9OnjxJUlISSUlJnD17li5duvDoo4/y+OOPExZWa9Vu+3D6kLpT//h28GqtbqQc+rSakr8RVJpM5CWtwPj++4iSEvwnTSLoqSfx6drVPna7KVJUbCA9Fffm0RUHMbqwqFiTdeESmzLOsunIWc6artCzQxuiRmh5UNe1yTdVAlRUVPDxxx+TmJjIjh07qKysZOTIkURFRTF58mRatWrlmIYv/ARfvgrfJIOoVIuI3fU0dLuzwZspQS0glpf4Dvlr1yKqqtAYHiRo3jy8m6K8gAsiRcUGUlTcG3cSFQuVVYJdx3KJ/yyLb84UENTOh0eH9mTWkB5o2jhn30tOTg4rVqzgnXfe4eTJkwQGBjJz5kwee+wxBgwY4JhGC3Pg63hIXwElBRA6GO56EvpOAk+vBl+2/LffyIuPJ399CoqioJk+jaDoaLwcvX/HxZCiYgMpKu7NoysOkl9cxhY3EhULQggOZBtJ+DyLT3+8QBsfT6aFd+Oxu3vRNcA5JXmrqqrYs2cPSUlJbNq0idLSUnQ6HXPnzuWhhx4iIKDWqt6No/SSWp3ywJtgzAb/bvCHGLVKZeuGp6EpP3uWC2+9RcGmzSg+PgTOmkng3Ll4OeIeXBApKjaQMRX3xp1FxZrjuYUkfJ7N1qM5CGDC7cFEjwijX4jji4fZwmg0snbtWpKSkjhy5AitWrViypQpzJ07l5EjR9p/Y2VVJfy0C756Qw3q+7RT4y5/iIHAhgffy06e5MIbb1K4fTsebdsS+OijBD76CJ7tGrcKzdWRomID6am4Ny1FVCzkmK6QtO8EHxz8leKySkbc3JF5I7TcFdbBqUkWjxw5QlJSEmvWrCE/P58ePXrw6KOPMmfOHHr06GH/BnOOqp7L9xtAVEGfcea4yx8aHHcp/flnLix/naLdu/H09yfw8ccInDkTjzbO8QqdjRQVG0hRcW9amqhYKLhczvtfn2LF/pNcvFTKbaH+PD68F/f164Kvj/PK8paUlLBlyxaSkpJITU0FYNSoUcydO5fJkyfbvyRyYQ4cTIDDK6DEBCE6db9LvwcaHHe58v0xLix/jeLPPsczKIig6Gg006fh4eOceJazkKJiAykq7k1LFRULJeVqjZeEz7PJvlhMWx9PRvfvwqSBoQwL64CXE3bsWzh16hTvvfceK1as4OTJk3Ts2JF58+bx5JNP0sXeK67KiuHo2t/jLgG91BxjA2Y0eKf+5YwMLix7lcsHD+IVHEzQE/PQTJ6M4qg9O80MKSo2kKLi3rR0UbFQVSU4cCKPrUdz2PndOQpLKghq58P420OYeEcIA7tpnDY9VlVVxSeffMLy5cvZtm0bXl5ezJgxg2eeeYaBA+2cn6uqSt3n8sV/1AzJ7bvCsAVqUN+7YUW+ig8c4MIry7jyzTd4d+9Ox6efov24cSiOTGfTDJCiYgMZqHdvpKhcS2lFJXt/vMCWo2dJ++E8ZRVVdA9sw6Q7Qph0Ryi9OzkvAP3zzz+zfPlykpKSKC4uZsSIETz77LNMmDDBvjnHhIBf9qji8utX0LajOi02+LEG1XYRQnDps8+48OprlP7wAz69w+g4fwF+EXq3TbkvRcUG0lNxb6SoXJ/CknJ2fZ/L1m9y2P/LRaoE3BrankkDQpkwIIQu/k2fEgbUgmJJSUm89tprnDp1il69erFgwQLmzp1L+/Z2XtF2cj98/iJkf6ouQf7DE+qKsTaB9b6UqKqiaHcqF5Yvpywri9a3306Xv/4Pvrffbl+bmwFSVGwgRcW9kaJSd84XlrD923NsOXqWb84UoCgwpFcHpoV34/5buzht1/6WLVt45ZVX2L9/P35+fjz22GPMnz8frVZr38bOpsPnL8GPO9TlyIPnqivG/DrX+1KispKCrdu48PLLVFy4gP/kyXT687N4dexoX5udiBQVG0hRcW8eSTqI6bIUlfpy4mIxW46qKWFO5V3G39ebyQNDmXFnd27p4ucUmw4dOsSrr75KcnIylZWVTJo0iWeffZbhw4fbNx702zH44mU4thE8fdS9LsP+BJpu9b5U5aVi8uLfJm/le3j4+BD05JMEzp6F4gYrxaSo2ECKinsjRaVxVFUJDmTn8cGh0+z6Ppeyyip03TVMv7M7428Ppo1Pw9OhNJScnBzefPNN3n77bf5/e3ceV1Wd/3H8ddhMREVMVFwySFxKUVxTQ0qcGtMWcy1LmxRzwy23Zvr96pFmYOIy7luUMlbSpDWTGS64TLnhElqZghtiiiIWiKzf3x/38BtCUC733AX4PB8PHgmc+z1fHh143+9+/fp1AgMDmTJlCgMHDsTNyD/W1xNhXyQc/8T0efuXoed0qOVjdlE5585x5f1wMuLicGvWjPpvzsIjKMi4utqBhEopJFQqt+HrDpKelcuWcd3tXZUKLy0zh38eSWbjwQskpmZSs5oLz7TzYWjnpjzSqPzboZTXrVu32LBhAwsWLODnn3/Gx8eHCRMmEBoaipeX+eMhpUq/aDrX5cjHoDmZtt7vMblcW+9n7NnDlffmknPuHB7BwdSfOQO3Zs2Mq6sNSaiUQkKlcpNQMZ5SisPnb7Dx4AX+/cNlsvMKaNOoNkM6N+GZAB9q3mfbdRoFBQVs27aNBQsWEBsbi7u7OyNGjGDixIn4+xt48NaNcxAXDj98Aq7u0HWsaet9M/cXUzk5pG2I5trSpRTk5FB3xHDqjn4dZ48axtXVBiRUSiFTiis3CRXrunkrl83HLrHx4AV+/vV33N2c6du2IYM7NSWwqe3XviQkJLBw4UI2bNhAbm4uffv2ZfLkyQQHBxtXl9RTsGsO/LjFdGhY94mm2WJu5oVCXmoqVyMXcPOLL3CpVw/vN6ZSq1+/CjMFWUKlFNJSqdwkVGxDKcXx5JtsPHCBr35I4VZOPr71ajCwQxP6Bzaivo1Pq7xy5QrLly9n2bJlpKam0q5dOyZPnsyQIUOMG3e5fBx2zobT30INbwh6AzqMABfzzpHJOn6cX+e8x+0ffqB6QAD1//ZXqrdpY0wdrUhCpRQSKpWbhIrtZWTn8XXCZTYdvsihczdw0qCnfz0GdmxCr1beVHOx3dTk27dvEx0dTWRkJD/++CMNGjRg/PjxjBs3Dk9PT2NucmE/7HjXtDNy7SbQc4Zp+xcz9hZTBQXc3PIlV+fPJ//aNWr16UO9yZNwa2L+jDNbqZKhommaL5CmlEov7RoJlcpNQsW+zl7LJCb+Ip/HX+LX325Tx92VZ9s1YmDHxjzsY7vBfaUUsbGxREZGsm3bNmrVqkVYWBiTJk2ibt26RtzAtHhyx7uQcgTqPgTBs+Dh/mBGd1Z+RgZp69Zx/cMoVF4edYYM4f6xYxzyDJcqFyqapoUAM4BwpdT20q6TUKncJFQcQ36BYu/pVDbFJxN78go5+QW0bliLgR0b81y7RtSpYbt1G8eOHWP27Nl8/vnneHh4MG7cOKZMmYK3t7flhSsFP//bNOZy9Ufwbm1qubR6xqxwyb16lWtLlpIeE4OTuzt1R47Ea/grOFUv3/5k1lDlQgVA07SVwCYJlapLQsXxpN/K4cvjKWw6nEzCpZu4OTsR0tqbYV0f4FFf2537cuLECebMmcOnn37Kfffdx5gxY3jjjTdo2LCh5YUX5MOJf8LucLh+GrwfhuAZ0LKfWeGSnZjI1cgFZOzYgYu3N/XCJlD7uefQXGy/Pqg4hwwVTdMGAJ2UUjNK+N50IAnwAlBKrSpH+RIqVZyEimP76fJvbDqczOZjl0jLzKFTszqE9WpOj4fut1m4nDp1ivfee4/o6GhcXV0ZNWoU06dPp3HjxpYXXpBvOihsdzhcPwP1HzG1XFr2NStcbsXHczViHlnHj+P2kB/eU6bi8biBM9rKobRQscvcNU3TQvTQGA3cMVqmaVo4kKSUitHDxE8PICFEJdKqYS3+p19rvpv5BO8++zDJN7J4ee1B+i//jl2nrmKLN70tWrTgo48+4tSpU7z00kssX74cPz8/xowZw/nz5y0r3MkZ2g6CcQfh+VWQmwWfvQwrg+Cnr0xb8ZeBe4cOPPDJRhotXgR5+SSPHcv5l18m6/hxy+pnBfZuqYQDnkqp0cW+fkMpVafI5yHADKVUb/3zUEoII+BI0VaJtFTEK+sOclNaKhVGdl4+MfHJLNuVyKX0LNo2rk3YE83p1crbZu/Kz507R3h4OOvWraOgoIDhw4cza9Ys/Pz8LC88Pw9OxMDuCEhLhAZtoOdM03HHZfz5VG4u6TExpC5ZSv7169R88km8J0+y+cp8R+3+uiNUNE0LBHYUC5VAIF4pZdZTJaEiXll3kN+yctksoVKh5OQV8M8jySyNO8PFtCwe9qlFWK/m9G5VHycn24RLcnIyERERrF69mtzcXIYNG8Zbb71lXLgkbII9EaaTKBu0Mc0Wa9GnzOGSn5FJ2ocfcv3DD1E5OdQZOtSmM8UcqvvrHryAtGJfSwfQNK3ME8v17rKOwEA9lIQQFYSbixNDOjdl59Rg5g1oS0Z2HqPXx9Nn8V6+TrhMQYH13ww3btyYxYsXc/bsWcLCwvj0009p0aIFr732GmfPnrWscGcXaDcUxh2C55ZDdgZ88qKpW+yXbWUrwqMG9SaM56Ft3+D5wgvciI4m8U9Pcn3tOgpyciyrnwUcMVQ80QfniygMmTLvEqePx3RQSo1WSh0xrHZCCJtxdXZiYMcm7JjSk8hBAeTkFTA2+ghPLdrDV8dTyLdBuDRo0IDIyEiSkpIYP3480dHR+Pv7ExoaavmYi7MLtHsRxh/Ww+V3+Mcg+MdgSCtbcLnUq0fDd97Gd8tmqrdvx9V580jq8zS/bd1qkzGp4hwxVEpaqFgYJsVbMBZLSUlB07T//5D9v4RwPC7OTvQPbEzslJ4sGtKOAgUTNh7lyYV72Jpw2SZ/PBs2bMjChQtJTEzk9ddf56OPPqJ58+aMGTOGixcvWlZ4YbiMOwi934Wze2FpF4h73zS4XwbVmjen6apVNFm7BqcaNbg0eQrnhwzl1pGjltXNTI46pvKH8ZPyjqmUhYypVG4yplI55Rcovk64zMLtv5CYmkmbRrWZ9mQLHmtuu6nIFy9eZO7cuaxZswZN0xg1ahSzZs2iUaNGlhd+8xJ8+zfTQWF1msGfI8D/yTK/XOXnc3PzZlIXLiIvNZWaTz2F95TJuDVtannddBVmTEXvqireWvECSh1st0RhS0VaKEJUHM5OGv0CfNg2KYh5A9qSlpnDK+sOMmTVfuLPG96hUaImTZqwbNkyTp8+zYgRI1i5ciV+fn5MnDiRy5cvW1Z47UYw8EN4ZYvp9Ml/DIKNQ03b75eB5uyM5wsv4PfNVu4fN46M3btJfLovV94PJ//mTcvqdq97O1pLpcjXDymlYkr63EjSUqncpKVSNWTn5fPJwYv8fecZrmVk80RLb6b+yd+m+4udPXuWOXPmEBUVhaurK6+//jozZsygQYMGlhWclwP7l5mmIat86DHFtN2+a9l3fs69cpXURYu4+cUXONeqxf3jxlJnyBCLjjV2qCnFendWCKbFj17AXGB70QH1IivqfYH08qyoLws5T6Vyk1CpWm7l5BH13TlWxCXy2+08+gX4MDmkOb71PGxWh8TERGbPns3HH3+Mq6srr776KtOmTcPX19eygm8mw7a/wo+boc6D0GceNO9tVhG3f/qJKxER3Pp+P64PNMVn7lzcA8s3OdahQsWRSEulcpNQqZpuZuWyek8S6/5zluy8AgZ2aExYr+b4eNpuQ8YzZ84wb948oqKiyMvLY/DgwcycOZO2bdtaVnDiLvh6mmlPsZZ94cn3oM4DZX65UorMPXu4GrmARgsiqVbOsJNQKYWESuUmoVK1pf6ezbK4M0TvvwDAsK4PMPZxP+73MO8gLUukpKSwYMECVqxYQUZGBk8//TSzZs2ie3cLnsm8HNi/VO8SK4DHpkK3CeBa9tBUSlk0qaHCDNTbmgzUV25V/U1TVVevZjX+t9/D7JoWzHPtfYj67iyPhe/ir18kcPrK7zapg4+PD/PmzePChQu8++67HDhwgB49evDYY4/x9ddfl+8ZdXGDHpNh/CHTrLBdc+DvHeDYxjLvJ2atWXLSUpGWSqX28toDZGTn8cVYaakISEzNYEVcIluOp5CTV0A3v7oM79aMkFb1cbbR9i+ZmZmsXbuWDz74gIsXLxIQEMDMmTMZMGAALuXd0v7cPtMU5JSj0KAt/Gk2+PY0tuLFSEulFNJSEaLq8KvnwbyBAeyf1YvpT7Xg3LVMRq+PJyhiFyt2J5J+y/rbm9SoUYOwsDASExOJiooiOzuboUOH0rJlS1auXMnt27fNL7RZDxi5E/qvgawb8PEzED0Irv5s/A9wD9JSkZZKpSYtFXE3efkFbP/pClHfnWN/UhrVXJx4rl0jhndrRmufWjapQ0FBAVu2bGHu3LkcOnSIJk2aMH/+fAYMGFC+Lqrc23BgBeydDzkZEDgcHn8TPAw42bIIaakIIUQxLs5OPPVIQz4JfZRvJj1G/8DGbDl+iT6L9zJoxfd8nXCZvPyyjVGUl5OTE88//zwHDhwgNjYWLy8vBg0aRK9evTh58qT5BbreBz0mQdgx6DQKjq6Hxe1h9zzIuWX8D1CMhIoQQgAtG9Ribv827J/Vizf7tCTlZhZjo4/wWMQulu46w2+3c616f03TCAkJIT4+nqVLl3Ls2DECAgKYNGkS6eklbYl4DzXqQp8IGHsAfINh12z4eyAcjTadSGklVT5UZExFCFGUp7sboUF+7J72OKtf6Wgah9l2iic+iOPTQxesvjOys7MzY8eO5fTp04waNYrFixfj7+/P2rVrKSjjzK4/uP8hGBINr26Fmg1hy1hY2dO03sUKqnyo+Pj4oJSSUBFC/IGzk0bv1vXZMLILX43vQbO6NZjxeQLPLt3HoXPW31+sbt26LF++nPj4ePz9/Rk5ciRdu3bl4MGD5SvwgW4wcge8sBayb8L65yDB8J2vJFSEEOJe2jSuzabXH2Xx0PZcz8hh4IrvmbDxKJfSy7YtvSXat2/P3r17Wb9+PcnJyXTp0oW//OUvXLlyxfzCnJygzQDT+S1/nmc6xthgVT5UpPtLCFEWmqbxTIAPO6b2ZGKv5nx78ld6zY9j4fZfyMqx3hhF4b2HDRvGqVOnmDZtGhs2bMDf35+FCxeSm1uOsR6XatAl1KwV+GVV5UNFur+EEOZwd3Nhcm9/dkztSUir+izcfppe8+P46niK1XdwqFmzJhERESQkJPDoo48yefJk2rdvz86dO616X3NU+VARQojyaFzHnSUvBvLZ6EepU8ONCRuPMmjl95y4ZN3zSgBatGjB1q1b2bx5M7du3aJXr16Ehoby+++22XrmbiRURKVnm803RFXV+UEvvhzfg/f7tyEpNZN+S/Yx8/MfSP0926r31TSNZ599lpMnTzJ9+nTWrFlDmzZt2LXLOrO6ykpCRQghLOTspDGkc1N2vhHMa90fJCY+mcc/iGPprjP8buX1LdWrVyc8PJx9+/bh5ubGE088wYQJE8jMzLTqfUsjoSKEEAapXd2Vv/VtzbbJQXT19WLetlP0CN/F4h2nrb54slu3bhw7doywsDCWLFlCu3bt2Ldvn1XvWZIqHyoy+0sIYTS/eh6sGd6JL8d3p1OzOkTG/kKP93eycPsv3MyyXri4u7uzaNEi4uLiyM/PJygoiKlTp5KVZf2pz4VkQ0nZULJSe3ntATKz8/inbCgp7OjEpZss2nGa2B+vUPM+F17t/iCvdX+Q2u6uVrtnRkYG06ZNY8WKFbRs2ZKoqCi6dOliWPmyoaQQQtjJI41qs/qVjvw7rAfd/e5n8Y7TdA/fyQfbTnEj0zrb7Xt4eLB8+XK+/fZbMjMz6datG2+++SbZ2dadQCChIoQQNvKwT21WvNyBrRMfI8j/fpbsOkOP8J1EfPMzaVYKl969e5OQkMCIESOYO3cunTp14ujRo1a5F0ioCCGEzbVqWItlL3Vg26QgHm/pzfLdifQI38ncrT9xPcP4lkTt2rVZu3Yt//rXv7h27RqdO3fmnXfeKd9q/HuQUBFCCDtp0aAmS14M5NtJQYS0qs+qPUl0D9/J21+eJPmG8WefPP3005w4cYLBgwfz9ttvs2nTJsPvIQP1MlBfqclAvahIzlzNYOXuRDYfu0SBgmcDfBjd048WDWoafq+4uDiCgoJwcipf20IG6kshU4ortyr+nklUMA95ezBvYAC7pz3OiG7N+Obkrzy5cA8jPzpE/Hljt9sPDg4ud6DcjYvhJVYwPj4+pKSk2LsaworKdc63EHbk41mdt/q2ZvzjD/Hx9+eJ+u4sLyz/ns7NvBgT7Edwi3oO+1xX+ZaKEEI4qjo13JgY0pz/zHyC/+3XmkvpWbwadYg/L9rLlmOXyMsvx0mQVlZpQ0XTtFD9Y6Wmab72ro8QQpSXu5tpwWTctGAiBwWQX6CY+Mkxgj+IY/3357ida93zXMxRKUNF07RA4LBSahWwSf8QQogKzdXZif6Bjdk2KYg1r3TEu2Y13tpykqCIXRy5cMPe1QMqaagAvsBo/d+H9c+FEKJScHLSCGldn8/HdOPT0K5Ud3NmyKr9fHnc/uPDdg0VTdMGaJoWXsr3puvfD9U0LdSccpVSMcAM/dMQYLuFVRVCCIejaRpdfOvyxdjutGvsSdjGoyzecdrqJ1DejV1CRdO0EE3TpmNqTXiW8P1wIEkpFaN3YflpmjbAnHsopdL1fw4GRllaZyGEcFReNdxYP7Iz/ds3IjL2F6Z8dpzsPPuMs9glVJRS25VSEcCRUi4J1VsbhWL5b3dW4SD89BI+QooWogfXqCIBI2xM1v8Ie6iKz101F2fmDwrgjT/588XRS7y0+oBVtny5F7uuqNdbJJ5KqaKBEQjsUErVKfa1eKVUmSdm6y2b7UqpdE3TQpRSJXaByYp669I0za5N8WFrDpCVm8/nY7rZrQ7C9uz93Nnbv35IYepnx6lf6z7WjejEQ94eht+jIq2o9wKKLx1NB9A07Y6uspLoIbQaiNc07Qb/HV8pEyPf5VhalrmvN+f6slx7r2sq8ztCW/9sRt9Pnr2Ky9KfrW9bHz4J7cqtnDyeX/Yf9p2+ZtX7FeWILZUBwOpiLRVP4Abgp5RKMrIOJbVUjHyXY2lZ5r7enOvLcu29rrH0+9ZmSUvF1nU3+n5V+dmz93NnKaPqn3zjFq9FHeZMagbvPvsIL3Zpatj9SmupOGKohACbioWKL5AI1DF6fETTtEzAvciXLuv/NWpuno+FZZn7enOuL8u197rG0u87MlvX3ej7VeVnryI/d1Axnr0HlFL1in/REff+SuPOGWGe8IcZXYZRStUwukxxd3prNB3T+qHtRrc+hbgb/U1qmkzgsQ6HG1NRSh1BH0MpwgtZa1Ip6L/QnfQZgKuAEtcpCWENek/ISuCObhthDIcLFd2qYutSemN6EETFNwBTV2ahQHtVRFQ9+ixQaRlbkV26v/TZWSGY/sB4aZqWiKkb5AiAUmpG4Yp6TF0kicXWrQgHoP//6aSUumN2nb5GKAlTKxO9VQJQl2K/1JqmeUpXhDBHOZ89YQN2CRU9PI4AEXe5ptTvCfvSuxACMbUg73jXp0/AOFT4RkDTtHBN0wYUeWPgZbPKikrFgGdPWJmjdn8JB2bhjgjXi13rJa0UUVaW7sYhrE9CRRhK79osLg1TdydADNBBv9YTmYAhDFKGZ0/YgCNOKRYV2113RFBKJWmaFl+kG8Os3Q6EuIt7PXvp+lhMR/1raYXjuMI4EirCaJ7cOWZS+IvuBaQXGTiVVoowUlmevRhMrWVhJdL9JYxW0vhI4S968XeRQhhJnj0HIKEijGbTHRGEKEKePQcgoSIMJTsiCHuRZ88xSKgIa5AdEYS9yLNnZ3bdpVhUTEV2RBiN6Z3gXIrsiKBfU7iq2Zc/Ds4LUW7y7Dk+CRUhhBCGke4vIYQQhpFQEUIIYRgJFSGEEIaRUBFCCGEYCRUhhBCGkVARQghhGAkVIYQQhpFQEcIGNE3z1DQtVtO0TZqm+dq5LtP1eky3Zz1E5SShIoTtHFFKDdTPlAnXNC1R07TEu71Av07p/zUkjPSTE0cBfkaUJ0RREipC2IFSagb6AWWlnFhYVJJSaoZS6o4z2YVwNBIqQthXDKWcoa6fjhlr2+oIYRkJFSHsayUwyN6VEMIocpywEGbSu6tWY9oJdyWm3XB7K6UGmluWPr6SpGnaAP2o28J7ePLfnXbvdm/0a/z0LrXidS3csRfAS3bsFdYmLRUhzKRvsz4D0x/zJOAzLOumWsmdXWAdSxpDKXLvQOCwUmq7HhSxmqb9oQ6apoVjGo+J0QPrs2JnjQhhOAkVIconDfBUSiUppSw6s0N/bYjeOinrvZOKHpGrlNoO+OrjMIUtndCirR8gFNOhVUJYjYSKEOVn5GysGEx/9Au7uA6Xo4wjmFowYDrIqnj9VqHPOBPCWiRUhHAMRbvAvIq2Qoyit6gML1eIoiRUhHAAevdVWdas3E0gsF3/d9FWixA2I6EihP0UXyEfA6wuDJh76Fh0DEYfgD9SeFa7Psi/StO00KIvkoF6YW0ypVgIM+mtiXBMg+vTgVXmdivpM7Oma5rmB8zQX7+y2DUDMHWJ+erXzy1yn8OYggXAE+hUfEqzUmq0vs9XKKbBfYoN3AthOE0pZe86CFHp6a2KWSWtJSlHWYFAuFKq3DO59PqEK6VKXM0vRHlJ95cQQgjDSKgIIYQwjISKEBVICeM5QjgUGVMRwgb0MYxNQDqmgXm7bWOvh1En4JB+tooQhpFQEUIIYRjp/hJCCGEYCRUhhBCGkVARQghhGAkVIYQQhpFQEUIIYRgJFSGEEIb5P0e1gGwhbkw6AAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x1185d9f50>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1, 1)\n",
"\n",
"__=ax.loglog()\n",
"__=ax.plot(rmids, xi_all_sats, label=r'${\\rm all\\ satellites}$', color=mred)\n",
"__=ax.plot(rmids, xi_noniso_cens, label=r'${\\rm non-isolated\\ centrals}$', color=morange)\n",
"__=ax.plot(rmids, xi_all_cens, label=r'${\\rm all\\ centrals}$', color='k')\n",
"__=ax.plot(rmids, xi_iso_cens, label=r'${\\rm isolated\\ centrals}$', color=mblue)\n",
"\n",
"xlabel = ax.set_xlabel(r'${\\rm r\\ [Mpc]}$')\n",
"ylabel = ax.set_ylabel(r'$\\xi(r)$')\n",
"title = ax.set_title(r'${10 < \\log M_{\\ast} < 10.25}$')\n",
"legend = ax.legend(loc='upper right')\n",
"ylim = ax.set_ylim(0.01, 2500)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Galaxy clustering signal shows differences in bias out to large scales. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compare halo concentrations"
]
},
{
"cell_type": "code",
"execution_count": 210,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAENCAYAAADgwHn9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3T2P49h6J/D/szAc7N5Arb2+QQPGDqi94SSqqszADZoKJ5OmP0FLkTGAAwmNTerCWBTUwQC7swk1n6Baym4ojbEXmKxK+gLtkmEYKBsDt5rBbGg8G/DwFItFipREvZD6/wChW+IheUiVzsPzRoqqgoiICAD+07EzQEREp4NBgYiILAYFqgwRaYrIUET6x84LPRERR0TcY+eD8mFQoK2YAnicsqwrIm3zOlgBraoLAHcAWofYnwlADyKiIlJbk84RkS8m7XBd2kMp+vsz63RFxDOv6DE2AYzNefoiIlMRaRZxHFS8vzh2BqhczI/5rXnrJCzvAoCqTsx7R0Q8Ve0dKIuLA+0HqjoQkc8AhgjORdq+XQA1AG9M4DqafXx/ItJV1VFsG3MAjfAzVX0lIjVV9Ys5EtoX1hRoI6q6UNUBgNuUJL1oAaGqSwSF4k7MVeh01+3swdK8XhSwACAibQArwNZk9iLv+Sn6+0uq9Zj16/EmIwaEcmBQoMKYAiKpWcDfpk1ZRGoi0jeF3UpVD9IslJc5pgWCoHCVsNyJLJvtYf+Fnp8tvz8HQLy5CFgTKOm0sfmIiuQASLoaXCEobHIVjKaAeW+256nqh10yZZozVpE8jsKr1si+7syypUnXytHk5ajqTETSCsCmqk5M2723yzFEFX1+Ijb+/lR1ISIXCbWA6LlELKg0EfkO6LQwKFCR6ngqfKN8AP81a2VzZT0w27kporlFRDwAQ9MMEhaoYzx1Rv+IoGCdmeUPqtoAMNlgNw8ALmP7dQHMzP4cFFBT2Mf5idnq+4vnwzSZLcNzCtPXEvkOlnj+HdAJYVCgozOF3RBBgWQL8AK22wRwGd2eqvoisox0jroA3sXXyypwTWEfFqDPagphU4rZV9v8f+sCfF/nZx8itZg34Wfx/Krq0nRgZ55nOjz2KVDR6gmf1QB8zrGuj+Qr1W1dItKEEfEA4ML8f4Xnea4juQklzsXT1f8CwTHaZZGr5BaK608o+vwk2eX7A4Lg1cnRNOQjVrui08CgQEW6x/PCMVTHmqGiqrpU1Q6CdvehGee+707KsPDzYEbXmH3e57wSr4cFX6RZJJykFQ0CLoCdRk0d8Pxs9f2FzJyGZzUZc06SbrC2wv4DHG2BQYEKYwrJZcJIlFrkynnd+kvTuTsA0BOR8Y4zYe+RPJqmgaeCeglgZZp53B1G8PjhvmKd2In9CaawbG+ygz2cn/j2t/7+TGf+JBYQXAQFf1KH/SUOOKeE8mNQoG0lNTMAQfPB+/CNadffqPlEVX0zlv4dgKaZAbtRAWq2swCwiM6eNQXeZWQs/pWqTsxrlLihmJTZuEsAb2OFpwvAT+iI7Uf+3970qr+g87PV92eC2TgaOEzhfx+pMdXCYJXUjGQCyMdT7hs5a6rKF1+5XwiufIcIZqwqgiaNbixNF0GB2AbQL2i/XQDTHHkbm3z1I5/3TV7a5v+1yDIXwBdzPHOzfnvNPjyT/iGazpwTJ5KPMF14jqL7bJr0YwS1k4OcnyK+v8j5ih6rpryix9w32+0X9TfB135eYr4worNjrtDbMGPmzdVvHUHzTHg1vo/9DhF03IbDPG+UY/bpRDAo0NkyhbOnCc0YIjLVPc2gjjQXuQA+ArwFBJ0O9inQObtDwn19TDv63jpBTRBaIWiH9xkQ6JSwpkBnzXSIRm/JUEMw3DRXpzNR1TAoEBGRxeYjIiKyGBSIiMhiUCAiIotBgYiIrNLfOvu3v/2tfvXVV8fOBhFRqczn839X1b+Kf176oPDVV1/h/v7+2NkgIioVEfnnpM/ZfERERFaumoK5q+MS5s6KeSb2mLs2XiXdPyZyl8grAHcaecasWc9B8DjEFYKbaE2SbkVARETFygwK5v4wd6o6Cd+LSDt8n5DeRXAXyBYSnnolIp5GHoguInMRQSQw1BHcxXGI4B717xgQiIgOI0/zUTcWAKZIfmgGAEBVZ6aAf3HvGHMXyvh9XjxE7t9uvALQUNVXacGHiIiKtzYopDxMZIWEm4jlVAfQT3ioyLMnPZmbhLF2QER0YFnNR3W8fI6qfdTgpnd3VNWliFzECvwXDzY3T2YKH6hei/Y5EBHR/mQFhfChI1FhkKjjZVNQJo08mtA0J7kALiJJZgBW+vScW09EurxrJRHR/mX1KSQV+mGQiNcgtjEG8CZac9Dg4eTR/U4RPAkr0ePjI0TEvq6vrwvIFhHRecqqKawQa+8P3+/6YBAzqmmYUHP4AuBVZPs+giGqiV6/fo3Hx8ddskJERMbamoIpsOOFfx2xPoBNmbkIU1WdmffRDu0PsYATfQAKUaG++eFnfPPDz8fOBtHJyDMkdWQK8VALwTBSAMHzZmPL1zLzGOoA7kWkZkYivQVs7eNzbJUO1jQfERFRcTKDgpmR7IhI28xEfojNHXARmbcgIk2Trg3gWxHphzUB0zw0RRBUvpjXA543D43MOt3Ig9U5V4Eqq9PpYDTafhzFbDZDo9FAr5c6fahQh95f0cqe/30r/eM4Ly8vlTfEo7hvfvgZ//jLr/b9f//dbxLThWnSlm/qT3/7NxuvM5vN4DgOHCe16yzTaDTCfD6H53nZiU36brd7sP0VZdd8R7dzjPyfEhGZq+pl/PPS3yWVqOxcd9u5oNvxfR8PDw8H3WcRyprvsuFdUuks/eMvvz6rSRyL7/tYLBZYLg8zlsL3fbx79+4g+ypSWfNdRqwpEB3RarXCYDCA4zi2KePDhw9oNpvwfR/T6dR+7vs+RqMRHMfBarWC4ziptYxooJlOpxgOh6jVapjNZjYQhfsJt7FcLuF5Hq6urnB3d4f379+jVqvZfd/c3ODq6goAcl2xR9ep1Wqo1+toNptr9zWbzTAYDOC6rt3XdDrF119/nZjv2WyGXq+HXq+HWq0Gz/Pw008/oVarpZ6DJGnn/BwxKNBZKLrvoCiO46DT6WA+nwMI2rrjBXXozZs3Nh0QdFBHC9qoTqeD8XiMZrOJ1WqFm5sbDIdDtNttLJdLfP78Gf1+/9k6rVYL8/kctVoNjuPg3bt3GI/Hdt9hYQsAd3d3mcf25s0bjMdjOI6DxWKBTqdjg0navlzXRa/Xg+d5GA6HAILA+fDwgFar9SLfYfrb21vM53PU6083YEg7B3Hrzvk5YvMR0QlxHAe9Xg+j0Qi+79tO1clk8qIj+u3bt7i5uUncznQ6tcHi8vISi8WLmxY/E9ZAwkK/2WxiNgumI4X/Rq+yG43G2u2F64R5bjabmE6nmfsKRY+1Xq/D99fPlQ3Tt9ttu9285yDtnJ8r1hSITojruhgOh/A8D71eD91uF57n4e7u7tlVMADbRJLEcRxMJhOsViv4vo/Vav1dacIr+Gjh/O233wIImqLi+86StE5YcK/bVyi+7mq1WhuIkkZu5T0Haef8XLGmQHRCZrMZ2u02ptMpvnz5guVyieVyiUaj8aJQ830/sekIAC4uLuA4Drrd7trRTZNJMAUobPd3Xde+woIxbH7ZxLp11u0rrzDf6+Q9B2nn/FwxKBBFHHtU0nQ6tVf/tVrNFvrdbvdFQXV7e4v37+PPp4LtTI526gJPI50cx3nRHBP2NUQ/D6/kwwI1uizat5HEdV34vv8sz+EEvXX7Wicp32myzkFU2jk/V2w+osoI72G0zQSyY1ksFvA8D77v25m20SvVRqNhm0bG4zEGgwGurq6wXC7R6/XQbDaxWCwwHo+xXC4xmUzQbrfRbDZt2334Go1G6Pf7aDabuL29tctD4/H42QijpGWtVssWzB8/fsTFxUVqG/x8PsdgMECr1QLwfD5G2r6i5yPsR/E8D/f393j79i1Wq9WzfM9mM9ze3sL3fTQaDZsX13VTz4Hrus/O17pzfo44o5kqIx4Ukm50l7Ys/nmZAgvRNjijmc4G73pKtD0GBTorDBhE6zEoEIHBgijEoECld4gCnX0NdC44JJVKh09LI9ofBgUiIrIYFIiIyGJQICIiix3NVBrH6Edg3wWdGwYFKi0W2ETFY1CgSvre/263DXi/Mdt5fnO8T38f/Pv7tIf19P68235PWKfTQavV2vp5A+FT0ra5K2oZ9le0Y+WffQp0Mk5hqOlf/8e/4K//41+OmodTFRZQ23JdF4PBYKN1wjurHmp/Rdkl36Fj5Z81BSLKZZeAsA3f93M9C/rUlDXfIdYUiChT+ByCQz18xvd9vHv37iD7KlJZ8x3FmgLREc1mMwwGA7iua58tMJ1OMRgM7D39fd+3zwVYrVZwHAeu6+Zatyir1cpuN2zf/vDhA5rNJnzfx3Q6tZ+n5TdJNNBMp1MMh0PUajX7kJzFYmH3E25juVzC8zxcXV3h7u4O79+/t89l9n3/2XMa8lyxR9ep1Wqo1+vPHs6TtK+0c//1118n5jvsH+j1eqjVavA8Dz/99JN9pGrSOUiSds6LxKBAJ+8Q/QzH6kdwXRe9Xg+e52E4HAIICuDo+zdv3jx70lmn00G9Xs+1blEcx0Gn07H5GI1GLwrqUFp+k55o1ul0MB6P7eM7b25uMBwO7dPZPn/+jH6//2ydVquF+XyOWq0Gx3Hw7t07jMdju++wsAWAu7u7zGN78+YNxuMxHMfBYrFAp9OxwSRtX2nn/uHhAa1W60W+w/S3t7eYz+fPnkGddg7i1p3zIrH5iOgERK/s6/W6fbpZ+PSxqLdv3+Lm5iZz3X1yHAe9Xg+j0Qi+79sRSXnyGzWdTm2wuLy8fPGozLiwBhIW+s1m0z7KM/w3epXdaDTWbi9cJ8xzs9nEdDrN3Fdo03Mfpm+323a7ec9B2jkvGmsKdHTHHnGU5NOBn9McvXIEYB96f3d392JZ2OSQtW6S0WiU+XzlPM1PrutiOBzC8zz0ej10u114npcrv1GO42AymWC1WsH3/bV5B56ag6KF87fffgsgaIqK7ztL0jrhsa/bVyjp3K8LREnnNe85SDvnRWNQoJPD21Q/aTQa9so1FH0g/aaKurqczWZot9tot9vwfR+dTgfL5XLj/F5cXODHH39Eu93GYrHA7e1tYrrw2dPh86mjfRTh/8NnT29i3Trr9pVXmO918p6DtHNedP8Rm4+IEpzKfIVut/ui7fj29hbv378/Uo4C0+nUXv3XajVb6G+S37AzOdqpCzyNdHIc50VzTNjXEP08vJIPC+zosqxakeu68H3/WZ7DOQbr9rVOUr7TZJ2DqLRzXrRcNQUR6QNYAqgDgKpmzswQkTaAK1V9Mfsia3vb7I/OTzhr+e9q/+vFsqTP8mwr73b+1CumFrNYLOB5Hnzft+3xnufh/v7eXmWOx2MMBgN75drr9dBsNnOtW5TovmazGRqNBpbLpS3EGo2GvWJdl9/xeIzlcmnz12w2bdt9+BqNRuj3+/YqPlweGo/Hz0YYJS1rtVq2YP748SMuLi5Sa0nz+RyDwQCtVgvA89pA2r7Wnfu3b99itVo9y/dsNsPt7S1830ej0bB5cV039Ry4rvvsfK0750USVV2fQGQI4E5VJ0nvE9K7AJoAWgCWqtrbZHub7u/y8lLv7+/zHi+doLQ+hbD5KG35uqCQV9btMFKDQixvbOqishGRuapexj/P03zUjRXIUwC9tMSqOlPVDwDShhFkbW+j/RGl+d7/zr6IKJ+1zUciktRotQKw1Xz3rO0VvT+iUDww7FK7IKqyrD6FOoJCOcoHABGpqeqmA6LXbm8P+6MSO8WhqmnYjERVkdV8FBbUUWGhvdmA4HzbK3p/RES0gaygkHRlHhbO62eZbLe9jff3+PgIEbGv6+vrLbJFlA/7KKjqspqPVgiu3qNqALBlU87a7YnIxvt7/fo1Hh8ft8gKUfHYjERlt7amoKoLvLx6rwPInsGxxfaK3h8REW0mz5DUkZmIFmoBsDfcEBEntnyn7eVYTmfqGE03afs8hafEEe1D5oxmVR2ISN8U1A6Ah9g8AhdAB0A42axpPmsDqIvIA4CZqQVkbi/H/ujMHaNNn/0IdC5y3ebCTEZLWzYCMIq8XyCYuLZundRleZYTEdF+8C6pdJZ2vfIv4hYbRKeId0klIiKLQYGIiCw2H1HpsROYqDisKRARkcWgQEREFoMCERFZ7FMg2gGf00BVw5oCERFZDApERGSx+YhOFoeaEh0eawp0MGl3Fq3ig2t4F1UqKwYFIiKyGBSIiMhinwLtHZtRiMqDNQUiIrJYU6C9YQ2BqHxYUyAqEJ/pTGXHoEAHxwKS6HQxKBARkcU+BTq4p+aV+VHzsU/ncIxUTawpEBGRxZoCHU3Yr/C9eR9eXfP200THw6BAJ6dq90EiKhMGBaI94igrKhsGBaI9elnrYccznTZ2NBMRkcWgQEREFoMC0QFxNjedOgYFIiKy2NFMR3PWQ0+9PwT/9v583HwQxeQKCiLSB7AEUAcAVR1tm15EPABDVV2mrNsG4ACYAFgB6AKYpKUnIqLiZDYficgQwFJVJ6Zwb5iCe9v0LoAHEdHYq2uW1wEMATwA+CezLQYEqpRPv/yKT7/8euxsEL2Qp0+hq6qTyPspgN4O6WcALgA0Iq8PsdrHKwANVX0V2xYREe3R2uYjEWkmfLxCcLW/cXoRqSHWdGRqCDfRFVTVB+CvzTkRERUuq0+hjqBQj/KBoIA3hfem6e06Jogs49sxgWJltldT1Q85joXo5J115zqVQlbzUQ2mszgiLPTjn2+Tvqeqs9hnMwAfY30S3YR1AQCPj48QEfu6vr5OS0p0urw/PI1IIjqirJpCUhNOWLjHawQbpRcRF0Fn8jMJncpTBB3PiSOeXr9+jcfHx6RFRES0oaygsEJw9R9VA2y7/y7pewBuox+YPocvAF5F0vsIhqhSSXDGLlF5rW0+UtUFXl791xE08eyavo1gLkPch1gAcVLSERFRwfIMSR3F5hm0AHjhGxFxYsvXpjfrhLWJZwHEBIPPsf13AAxy5JOodHgvJDo1mTOaVXUgIv3ITOOH2NwBF0HBPcmZPrREcr/EyMyI9hHMYfA4V6Gc+HhNovLJdZuLdUNCzQihUeyztUNITY2gsWYZh6CWEK94icqPd0klIiKLQYH27nv/O07aIioJBgUiIrIYFIiIyGJQoO3x1gzF4zmlI2NQIDp1DBR0QHwcJ20tfEjM73OmZ2fzS/ac/O43x80IkcGgQHSKWDOgI2HzERERWawpUOHYTERUXqwpEJ2QT7/8avtqiI6BQYGIiCwGBaKy4NBUOgD2KdDGwruhfh97T0Tlx5oCERFZDApERGSx+YjoBHDEEZ0KBgXaGeclEFUHgwLRCYrWHH7P+yLRAbFPgYiILAYFIiKyGBSIiMhinwJRSWz6/AqibbCmQJm++eFnzlomOhOsKVBuDAzH8WIOQ3j/o96fD58ZqjwGBdoY5yUQVRebj4jKindNpT1gUCAiIotBgYiIrFx9CiLSB7AEUAcAVR1tm15E2gAcABMAKwBdABNVXW67PyIiKkZmTUFEhgCWqjoxhXPDFOzbpq8DGAJ4APBPJu1yg/WJiGhP8tQUuqo6iLyfAhgguNLfNv0rAPVoMNhhf7QHHH5aQhyqSgVYGxREpJnw8QqAu0t6VfUB+Lvuj4iIipVVU6gjKJSjfAAQkZop3DdOLyJdk64OoKaqH7bcH9HZiU9m420vqEhZfQo1mM7eiLDQjn+eN/0MwMdYn0F3y/3h8fERImJf19fXqQdD2/ne/44T1ojORFZNIenKPCyc41f0udIn9CNMEXQ8j7bYH16/fo3Hx8ekRUREtKGsoLBCcPUeVQNsv8BG6UWkBuALgFeR9X0EQ1S32R8dEGsLRNW3Niio6kJEkvoNZjuk/xAr4B0EcxI23h8Rgbe6oELlmdE8is0TaAHwwjci4sSWp6Y3weBzbPsdBENOc+2PiIj2J3OegqoORKQfmYn8oKrROQMugoJ9kjP9yMxY9gE0AHjR5TnWJ6I1wjkmf/rbvzlyTqiMct3mIjJkNGnZCEEncd70PoDU5Vnr035x0hrReePzFAgAgwERBRgUiCrmaZTY/Kj5oHLirbOJSu7TL7++fGQn0ZZYUzhzbDYioijWFIiIyGJNgZ7hrGWi88agQAAYDCqNz1mgDbD5iIiILNYUzhQ7mKvPzmz+yyNnhEqFQYGoIuLDUm2T4O9+c4TcUFmx+Yjo3Hh/4J1VKRWDQsV988PPbCoiotwYFIjOFWsMlIBBgYiILAYFonPHGgNFcPRRxWz6gBVOWiOiKAaFM8GncRFrA5QHg8KZCYMDawhElIR9CkREZDEoEBGRxaBwJr73v2OT0Zkq/MlsHK1UaQwKRGciMziwsCewo7ky0m5lEa8dsLZAmfj8hbPGmgIREVkMCkSUjM1JZ4nNR0SUDwPEWWBQqCg7Se3I+SCicmFQqCh2KFOWcCTS7/lkNopgnwIREVkMCkREZOVqPhKRPoAlgDoAqOpol/RmOQBcAbhT1Q+RZW0ADoAJgBWALoCJqi7z5JWIiLaXGRREZIig4J6E70WkHb7fNL2IeKrai6SfiwgigaEOYGhePoB3DAhr2BEh/xMA+xKIaDd5mo+6sQAwBdBLS7wuvYjUEBT0UR6A97HPXgFoqOqrtOBz9jiGnLYUv91F4fdGolJbW1MQkWbCxysA7pbp6wD6prYQvfqvRVdQVR8vgwetwRoC7U3axQdvh1FJWc1HdQSFepQPBFf9pvDeJP1SRC5iAaEFYBZdQUS6Zjt1ALVonwMREe1PVvNRDaazOCIs9OOf50qvqotwgWlOcvG8OWoG4KOqTkwHdcMEiUSPj48QEfu6vr7OOKRqYdWfisK/JQKyawpJTThhoR+vEWyTfgzgTbTmkNCpPEXQ6Zw44un169d4fHxMWlRp/PHSyWFzUiVk1RRWiLX3h+8Tmo42Sm9GKQ3jNQcRUVODCPkIhqgSEdGera0pqOpCRJL6DWa7pDdzEaaqOjPvm5Hg8CEWQBwEcx6I6AA2vv0FR8FVSp4hqSNTiIdaCIaRAgBExIktz0rvIggU96Zm4AB4C9jaxOfY/jsABnkOptI4BJXKgn+rpZY5eU1VByLSj8w0fojNHXARFNyTrPSmWWhq1vMi24hub2RmPPsAGgA8zlUgOjzeMO885brNxbohoWaE0Cj2WWJ6UxOQjH35ADgENcb+QI+cDyKqNt46u2xYLaeyiP6tckRSaTAolAyHotKxsVmp2njrbCIishgUiGj/OCKpNBgUThV/RHQi9nr7C/6dnxz2KRDR4fCOqyePQYGIcuEgh/PA5iMiIrIYFIiIyGLz0Yn55oefAQB/+ssjZ4TomNL6GNj3sHcMCicmfKzmp/CDv784Wl6I8ih0Mlu8I5ojkw6OQeHQeKVDFXESHc/8PRWOfQpEVH6bznfg/IhUDApEVB0s7HfG5qNNsKpKlFu0eYk3zysPBoUjeRpl9D+CDxhoqOR499RqYFA4knCUEfgDoopJ6oCOB4ydAwibiPaGfQpERGSxpnAiPnE+AlFx4v1/Wf2B7C+0GBQOxfzRncTYbqJzlXdyXJ7mqW0DyIkHIAaFA0kLBgwSRHRKGBSKlFZlJSIA2RdBhY9g2udvML7tvE1VJ45BYQP2D3bD9ES03sGDRZp9BJGSBQkGhT0I5yB8f+R8EJ26vMGgEkrScsCgQERUhE0L/W1vD77nmgeDwi5Svhw7MY2I9qKSs6dPpCbBoEBEJ6dSzUZZsoLBgfskGBSKwDkIRAdxMh3Seex65X+kmgNvc0FERBZrCjtgzYDouEpVc9jVgWoOuYKCiPQBLAHUAUBVR7uk33X50Z1IhxAR5ZMVHNYtr1RgySGz+UhEhgCWqjoxhXNDRNrbpt91eZGur6+3Wu/TL7+Wtpbwv//h8dhZODgec/XFjzftNxp+XubfcGjb8itLnj6FrqpOIu+nAHo7pN91eWH++Mc/7mOzJ+3//N9/O3YWDo7HXH3bHm88OMSDRtKyvO833fem9lV+rW0+EpFmwscrAO426XddTkR0bPGCPOt9KP6AoXj6ePNUPN2hmq+y+hTqCArlKB8ARKSmqv4m6XddnrC/g+CzDohoV7ve0uNQfRuiqukLg7b8H1X1VeSzGoAvABqqutwkPYDmLsvj+zPL/x+A/xz56F8B5G1Qfb1B2qrgMZ+HczvmczteYPdj/m+q+lfxD7NqCklX5nXzb/yKPk/6XZe/oKr/JelzIiLaXFZH8wpALfZZDQBSmnKy0u+6nIiI9mhtUFDVBV5evdcBzLZJv+tyIiLarzxDUkexeQItAF74RkSc2PK16QtYTpRIRNpmnkvSsr5Z3hWR7qHzti9px2w+75vfZy38/zHySOWytqPZJnqaYewA8KMzjM0PrKOqrTzpi1i+q5OfMV0wE2QdABMETXRdAJOkjvsyEhEXwSCFFoKJj73Y8iGAu3D+S/x9GeU45i6eLqZ8AO/KfLxR5vcLAFcIvscPCcsr9fted8yF/75V9axeAIYA2mnvq/gyfyRqXl+qerzmu/QSPv8Se+8CmB47v3s+5i6C/jjn2Hks+Hi92Ps5gH7sfFTq953jmAv9fZ/jXVIPNmP6xLxCMKz3lVbkijGPc54Qqaq+VqQ2CNjh6fE+Rw/A+8j7Sv2+cx4zUODv+6yCAguI6hQQG8iaMFlZpv8k7EfpZ69x8uoAkvpGakBlf99rjzlU5O/73G6dfZIzpg/BtDGvEJyDmsbaYSssnCkfFf4N1JE8N6YKZgBW4d+0iHgi0tUSt6+r6lJELmKFXwtPoxMr9/vOccwAiv19n1tQYAGBahQQG9h4QmQVJFw1ThG0r5f6O9dg2DoAW9NzAYT3oank7zvjmIGCf99n1XyEMy4gYldJUwCDY+XnwM5uQqQZgqqx5jEfwQiVKhkDeBMJgOfw+44fc+G/73MLCiwgAlUsIBLp+U6I/BD7m3YQDNOsBDOseBi9ikbFf99Jx7yP3/dZBQUWEFalCogczmpCpPmuP8c+7qAitUPXI0bNAAAAzklEQVTzXU5VdWbeN4Fq/77Tjtko9Pd9bn0KgCkgIsO2Kl9AiEhlCwjA/kBcAG0AdRF5ADALr6hUdRDOaEbwg3ko+7DcrGNG8HfeR1BINhCMdS/1MQN20l4dwCxyu/23AKLHXanf97pj3sfvO9eM5qrZ94zpU2P+kLp4KiBKPZuXzlPkNvpxE1XtRNJV5ved55iL/n2fZVAgIqJkZ9WnQERE6zEoEBGRxaBAREQWgwIREVkMCkREZDEoEBGRxaBAREQWgwIREVkMCkREZP1/b0sDpDaMykQAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x119f18d50>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1, 1)\n",
"\n",
"conc_bins = np.linspace(0, 25, 150)\n",
"__=ax.hist(iso_cens['halo_nfw_conc'], bins=conc_bins, normed=True, \n",
" alpha=0.8, label=r'${\\rm isolated\\ centrals}$')\n",
"__=ax.hist(noniso_cens['halo_nfw_conc'], bins=conc_bins, normed=True, \n",
" alpha=0.8, label=r'${\\rm non-isolated\\ centrals}$')\n",
"title = ax.set_title(r'${10 < \\log M_{\\ast} < 10.25}$')\n",
"\n",
"legend = ax.legend()\n",
"# xlim = ax.set_xlim(0, 25)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These distributions don't look especially different, but the long tail at high-concentration is a highly biased population, so upweighting by that tail gives a significant boost to clustering, even on large scales. \n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compare the galaxy lensing signals"
]
},
{
"cell_type": "code",
"execution_count": 179,
"metadata": {},
"outputs": [],
"source": [
"from halotools.mock_observables import delta_sigma\n",
"particles = return_xyz_formatted_array(\n",
" halocat.ptcl_table['x'], halocat.ptcl_table['y'], halocat.ptcl_table['z'])\n",
"particle_masses = halocat.particle_mass\n",
"total_num_ptcl_in_snapshot = halocat.num_ptcl_per_dim**3\n",
"downsampling_factor = total_num_ptcl_in_snapshot/float(len(particles))\n",
"\n",
"rp_bins = np.logspace(-1, 1.5, 20)\n",
"period = model.mock.Lbox\n",
"rp_mids, ds_all_cens = delta_sigma(pos_all_cens, particles, particle_masses, \n",
" downsampling_factor, rp_bins, period)\n",
"\n",
"rp_mids, ds_iso_cens = delta_sigma(pos_iso_cens, particles, particle_masses, \n",
" downsampling_factor, rp_bins, period)\n",
"\n",
"rp_mids, ds_noniso_cens = delta_sigma(pos_noniso_cens, particles, particle_masses, \n",
" downsampling_factor, rp_bins, period)\n",
"\n",
"rp_mids, ds_all_sats = delta_sigma(pos_all_sats, particles, particle_masses, \n",
" downsampling_factor, rp_bins, period)\n"
]
},
{
"cell_type": "code",
"execution_count": 207,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEkCAYAAADjOHzWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsnXl8VOX1/99PdkJCJossYZ8AKopKFgQ3VBKXWrVqAmprrd8WotX6bX9tCfhti3YRg7a1336rTbB1q2IIKlq1aoIL7kICbqgsA4gsAplMNrLOPL8/7p1hEib7TCYTzvv1uq+5c7fnzA3czz3PeZ5zlNYaQRAEQfAnYcE2QBAEQRh6iLgIgiAIfkfERRB8oJRKV0oVKqUWB9sW4ShKKatSKjvYdgjdI+IiBB3zQV7ayb5FSqlccxmwB73WuhLYAOQMRHumkO1QSmmllKWL46xKqWrz2MKujh0o/P33M89ZpJQqMhfv35gOlJr3qVopVaaUSvfH7xD8S0SwDRCOX8yHwgLzq9XH/kUAWus15nerUqpIa50/QCZWDlA7aK0LlFJVQCHGveis7WzAAswzBTBoBOLvp5RapLUu7nCNCiDNvU1rnaiUsmitHf75JUIgEM9FCBpa60qtdQFQ0skh+d4PGq21DePh2i/Mt+Ky/l4nANjM5ZgHNYBSKhewg8ezCgg9vT/+/vv58sLM85M6doWJsAx+RFyEQYn5oPHV3eHoS5+7UsqilFpsPjTtWusB6e7qKeZvqsQQlywf+61e+8oD0L5f708f/35WoGM3GHQhuMLgRbrFhMGKFfD1dmrHeGj16AFrPqiWmtcr0lqv6I9RZjeN3cvGYvdbtFdbG8x9NvO4nB505Vm11uVKqc4epOla6zVmbKOoP7/BG3/fHy96/ffTWlcqpTJ8eCXe95IO4pSO199AGDyIuAiDlSSOPsS9cQDJ3Z1svukXmNdZ7o9uJKVUEVBodu+4H8ylHA36r8R4QJeb+3dordOANb1oZgeQ2aHdbKDcbM+KHzyXQNyfDvTp79fRDrMr0Oa+p5ixKK+/gY32fwNhkCDiIgwpzIdmIcaDzSMEfrhuOpDpfT2ttUMpZfMKQmcDCzue192D2xQN94O4nefi7iIy28o11/ssBIG6P4HAy6ua597W0V6ttc0cKNDtfRYGFom5CIOZJB/bLEBVD8514PvNua9k4tU148UOIMNct9Pe5iR8dw11JJuj3kglxm/07PN6a8/Bf/EWf98fX/Tn7weGCOb1oMvLQQdvTwg+Ii7CYGUj7R+ybpLoYoiw1tqmtc7DiEsUmvMkAh0Mdj9EizBHQ5ltbuyhZ5DkfoB6dfe4Jwt6i0k20K9RbgN4f/r093Njzolp51mZ98RXMkQ7gRdKoZeIuAiDEvNha/Mxcsji9Sbf1fk2M4heAOQrpUr7ObN7I75HP6Vx9IFvA+xm91V2P0ZcOdxtdRgs4DPeYj50c3vTQADuT8fr9/nvZw6aWNNBWLIxBMTXwIhMBnBOktAzRFyEwYCv7hMwukWWur+YcY9edQtprR3mXIyFQLo5o7tXD2LzOpVApfdscPPBmek1lyNLa73GXIp9XqgDncwutwELOjyEswGHj4D3Yq/13N56IX66P336+5miWOotQKaIbPTy4Cxu0fPVPWYK0erBHDs6XlGScl8IFuaDMB/jwZkOFAMVPmZo2zC6WKz+GCprXjOvK8/CK/CdCxS42zUf5u4HWcehyNkYI5fc+21AiXuGuo82ioD5GG/kBV4z2QsxRp3ZvEZ1zce4B8Xmse423bPk3UOJ/TGSrNv7Yx7Xr7+f1/3K8PqtOzppLtHrNy/G8O7cgx38NXxa8CMiLoLgB8wHYy6m2Jhv40kYwuD2DgLRbiFGgNw9vHe5zPkQBgMiLoLgB7y9DR/7ygKVEcCrGywbWA2SGkUYHEjMRRD8wwZ85M0yu60CFmw2xcyOEadwiLAIgwXxXATBT5gxBO9UJRaMYcY9Cu4LwlBCxEUQBEHwO9ItJgiCIPgdERdBEATB74i4CIIgCH5HxEUQBEHwOyIugiAIgt+Rei4mKSkpetKkScE2QxAEIaSoqKg4rLU+oeN2EReTSZMmsXHjxmCbIQiCEFIopXb72i7dYoIgCILfEXERBEEQ/I6Ii8m+fftQSnHnnXcG2xRBEISQR2IuJqmpqezbty/YZgiCIAwJxHMRBEEQ/I6IiyAIguB3RFwEQRAEvyPiIgiCIPgdERdBEATB74i4CIIgCH5HxEUQhjjl5eWkpaWRn5/v83sw6YlteXl5FBdLpehQQ8RFEIY42dnZFBQUdPrd3/RGCHpiW35+PtnZ2X1uQwgOMolSEAS/4XA42LFjh1+v2VFYAtGG4H/EczGR9C+C0D8cDgcLFy70+zUrKyux2WwBa0MIDOK5mEj6FyHU8X4Il5WVUVhYiMVi6fP1VqxYQXp6Og6Hg7KyMoqKirpsp7y83CMG7nPdXofNZqOoqIisrCw2bNjA0qVLe2Sb3W6noKAAq9VKUVFRn9vo7LcIgUPERRA64cDdd9P8+RdBaTv65JMYfccdvTonLy+P0tJS0tPTsdvtLF++nMLCwj61X1xcfMyDu7t2cnNzsdlsVFVVsXjx4nbXy8nJoaKiAovFgtVqZeHChZSWlnZrh9VqJS8vj4qKCoA+tdHVbxECh3SLCcIQoaysjPT0dAAyMzOprKzs87WsViv5+fkUFxfjcDhYtGhRn9spLi7GarV6vIj09HTKy8v7bFtv2+jqtwiBQzwXQeiE3noOwcZqtbJmzRrsdjsOhwO73d7na2VnZ1NYWEhRURH5+fksWrTI05XU23bcwXdvQZk/f36fbettG139FiFwiOciCEOEjIwMrFYrixYtOmaEVW8pLy8nNzeXsrIyqqursdlsnu6knrazZs0aALKysrBYLGRnZ3sWfz3ce9JGV79FCBwiLoIwBHAHut3dVd6jq/rSPVZWVuY5z2KxeK7bXTtWqxWHw9HuWu44iff2/nSL9baNzn6LEFjCZeitQXFx8Z3SFyuEKlarlbfffhuHw0F1dTWpqal89tlnHDhwAKvVyl/+8hcqKioYO3YsTU1N7b5Pnz79mOvt3LmTpqYmtmzZwpYtWxg2bBjZ2dldtpObm8v06dN57LHHcDgcjBkzBqvVChjB9vvvv5+amhq2bNlCamoqY8aMobKyskvbmpqaWLZsGVu2bOGUU07BarX2uo3OfovgH+666679d9555zGzWpXWOhj2DDoyMzP1xo0bg22GIAhCSKGUqtBaZ3bcPuQD+koptzuSARRqraWzVRAEIcAMaXFRSqUDG7XWlUqpbKAUQ2QEQRCEADLUA/pWwJ1edaP5XRAEQQgwIeO5KKVygSyt9THpXJVSiwEbkASgtS42P9copdzDUrIB/87cEgRBEHwy6D0XpVS2KR75wDHJiJRShYBNa73GFJU0U4gA0Fq7xyYuACTjnSAIwgAw6MVFa12utV4BdDZYf5HWeo3X9zKOdoUBHs9moZfQCIIgCAFk0ItLV5gB+47YMbrA3MfkAsVaa4cZ1BcEQRACTMjEXDohCUNMvHEAKKUsGAH8lYBdKZWEEdSXuIsgCEKACXVxsWAG8b1wi02S1roSSBxYkwRBEISQ7hbD9FI64BabXqWEdVeidC+SFkcQBKHvhLq42Dl2BJkF2o0S6xGpqalorT2LiIsgBIa8vDyKi49JRdVjysvLSUtLIz8/v/uD/cBAt+dvgmV/SIuL2e3VUUSSkLiKIAxa8vPz+5U4Mjs7m4KCY6a7dUl/xKwv7fmL/tjtJlj2h7S4mBR7z2sBcoBeF4twd4uJxyIIgcWdXXmgcDgcnmJioUSo2u1m0Af0zeHG2UAukKSU2gGUm14LWusCpdRiU2CswI4O8156RGpqKvv27fOn6YIgdMDhcGCz2Tx17geivYULFw6omPmDULXbm0EvLqaIVAIrujim032C0Gf+swQOfBKctkfPgEvv6fHh5eXlFBQUkJ2dTVZWFmAUySooKPA8oBwOh6fWvN1ux2q1kp2d3aNz/YXdbvdc110pcsWKFaSnp+NwOCgrK/Ns78xeX1RWVnoKl5WVlVFYWIjFYvEUN6usrPS0476GzWajqKiIrKwsNmzYwNKlS7FYLJ62ly9f7rkfPfEgvM+xWCwkJSW1K6rmq63O7v2MGTN82l1eXk5+fj75+flYLBaKiopYt24dFoul03vgi87uuT8Z9OIyULi7xZYtWyZdY0LIkZ2dTX5+PkVFRRQWFgLGg9z7+7x586ioqPCck5eXR1JSUo/O9RdWq5W8vDyPHcXFxcc88N10Zq+vSpJ5eXmUlpaSnp6O3W5n+fLlFBYWeipUVlVVsXjx4nbn5OTkUFFR4fGiFi5cSGlpqadt90MbYMOGDd3+tnnz5lFaWorVaqWyspK8vDyPKHXWVmf3fseOHeTk5Bxjt/v4kpISKioqSEo6OhOjs3vQka7uuT8RcTGRbjHhGHrhOQwWvD2NpKQkz8N5zZo1x3ghCxYsYPny5Z4HamfnBtre/Px8CgoKmD9/Pu5qsD2x15uysjLP8ZmZmT6P8cbtEbnFIz093VMW2f3p/daflpbW5f1wn+O2IT09nbKysm7b8r4Pbtz3Pjk5udP23Mfn5h4NN/f0HnR2z/2NiIsgDCG832TBeAsG48274z53V0p35/qiuLi4W/HpSbdadnY2hYWFFBUVkZ+fz6JFiygqKuqRvd5YrVbWrFmD3W7H4XB0aTsc7ebyfsjPnz8fMLrYOrbdHb7Ocf/2rtpy4+vep6Wlddqer/va03vQ2T33NyIugnAckJaW5nmTduNwOHx2MfUEf73tlpeXk5ubS25uLg6Hg7y8PGw2W6/tzcjIYOXKleTm5lJZWUlJSYnP49asWUNubi5ZWVnYbLZ2MRz3enp6eqfnd0ZX53TVVk9x290VPb0Hnd1zf8fXuhyKrJT6hVLqwX4uv/CrxQFChiILQ5lFixYd07deUlLC0qVLg2SRQVlZmccbsVgsHvHojb3uoL138BzwBMStVisOR/vpcO5YjPd2t2fhfvB77+vOS8vOzvaMhHPjnqPSVVtd4cvuzujuHnjT2T33N915LrPofw2U/s8CGgAk5iKEMpWVlRQVFeFwODzxiqKiIjZu3Oh56y0tLaWgoMDzJp2fn096enqPzg2Ene6Z4zabzfMwTEtL87xBd2VvaWkpNpvNY196erontuFeiouLWbx4scercO93U1pa2m5EmK99OTk5ngf86tWrycjI6NRrq6iooKCggJycHKC9d9JZW13d+wULFmC329vZXV5eTklJCQ6Hg7S0NI8t2dnZnd6D7Ozsdverq3vuT5TWuvOdSv1Sa31vvxrwwzUGgszMTL1x48ZgmyEIghBSKKUqtNaZHbd32S3mD1EIBWERBEEQ/ItPcVFKnWHGWy702paglLpGKXXGwJk3cEjMRRAEwX8cE3NRSi3EKBO8EbhWKVWltb5Ya10DPK2UcgLhA2xnwJGYiyAIgv/wFdBP8+4/U0pZlVIPAgVa61pADZh1giAIQkjiq1usXRIdrbVNa30LsEgplQB0PgJAEARBEPAtLnal1DylVIlSaoR7o9b6PiCTIeq5SMxFEATBfxzTLaa1flopNRlYbXaDee9bp5TKGDDrBhCJuQiCIPgPn6PFtNY7tdZPd7Jvkzlq7EFvz0YpdXWgjBQEoe90LHMb6mV7YWj8hqFOrytRKqV+hFFGuF1uBq31MyIwgjD46FjmNphleyG0S/cKPacviSurzWHJviZHDsl4jCAI/iHUS/cKPafXngtdjxZL7KshgiAMbdyle4Xjg76IS5qv7i+zu8x3Tc0QQEaLCaFOZWUla9asYc2aNeTn5/c4o25nOBwOCgoKWLNmDeXl5e2y69psNs++goICT1vl5eVkZGR49rltsdlsx5Qc9i7OlZaWxooVKyguLiYjI8Nzvd78Jvc13ccKwaXX3WJa63uVUvcopR7i6JyYdKDYnA8TkshoMSHU6WmZ257iz7K97u++Sg6HUuleoef0qViY1nqJUmo54M4pXam13uk/swQh+Pz0pz9l8+bNQWn7jDPO4P777+/VOb0t9dsVgSjb2x2hULpX6Dn9qUSZAUw217f7wRZBEPpBb0v9dkUgyvZ2RyiU7hV6Tp/ERSn1KmAF3J2wN5s5/Rf4zTJBCDK99RyCTU/L3PaEQJfthdAs3Sv0nL7Mc1kOFGqtp2it55vLFGC1uU8QhAGmN2Vue0IgyvZC6JfuFXpOXzwXm9Z6XceNZtoYP5gkCEJv6U2ZW6vVekyZYF/4u2yvu62OJYdDqXSv0HO6LHPs8wSl5vkSl+72DXakzLEgCELv6VOZ407Q3jnFvBoYQYcJlkqpX/Th+kFB5rkIgiD4j750i90MzFRKdRxInglsVEq5E/4oYB5wXz/sGzBknosgCIL/6Iu4WDEEpruxhQq4pw/XFwRBEEKcvojLQq31pp4c6OXFCIIgCMcRvY659FRYenusIAiCMHToUlyUUn9XSv3IVwBfEARBEDqjS3HRWt8MVAMPKaVKpBiYIAiC0BO6jbmY5Y6fBlBKLVRKrQaqgFKt9WsBtk8QBEEIQXoVc9Far9RazweWABlKqVeVUsuVUmcExjxBEAQhFOnLJEq01jVa63u11hcBxcC1SqlXlFK/UEpN8qeBgiD0jLy8vH7Vp3cX7RqoQlsD3Z6/CXX7A02fxMUbrfVOrfUSrfXFwDpgiSk0MhBAEAaQ/Pz8PmUndpOdnU1BQe9mD/RHzPrSnr/oj91ugml/KNBvcfFGa71Ja32zKTQ7gRVKqQf92UagkPQvQqiTnZ09oMkaHQ6Hp7ZLKBGqdoca/SkW1iVmAsuQSWIp6V+EUMadHt9ddngg2lu4cGHIZR4OVbtDkYCJiyAIA4fdbqegoMCT4h5gxYoVpKen43A4KCsr82x3OByeNPZ2ux2r1dppd1plZaUnjX1ZWRmFhYVYLBZPrZXKykpPO97164uKisjKymLDhg0sXbrUUxLZ4XC0S9XfEw/C+xyLxUJSUlK7Gi++2iovL6egoIDs7GxPW2VlZcyYMcOn3eXl5eTn55Ofn4/FYqGoqIh169ZhsVg6vQe+6OyeH4+IuAjCEMBqtZKXl+epVV9cXHzMA9/NvHnz2tW0z8vLa/fA9iYvL4/S0lLS09Ox2+0sX76cwsJCT8GwqqoqFi9e3O6cnJwcKioqPF7UwoULKS0t9bTtfmgDbNiwodvfNm/ePEpLS7FarVRWVpKXl+cRpc7ays7OJj8/n6KiIgoLCwFDgHfs2EFOTs4xdruPLykpoaKiol2Z5s7uQUe6uufHI/0WF6XUZCAdI92+Mj/Ltda1/b22IASTu/79GVv2Beef8fTUESy7/JQ+n2+1WsnPz6egoID58+d7CnC5C3h5s2DBApYvX+4RAG/Kyso8x2dmZvo8xhu3R+QWj/T0dE+VSven91t/WlpaO6HriPsctw3p6emUlZV125b3fXCTlJRERUUFycnJnbbnPt67gFpP70Fn9/x4pV/i4p7fYk609N5+jVKqTARGEIJDdnY2hYWFFBUVkZ+fz6JFiygqKmLDhg3t3soBT9ePL6xWK2vWrMFut+NwOLDbu06G7vYovB/y8+fPB4wuto5td4evc9wP+q7actPxXLvdTlpaWqft+YrF9PQedHbPj1f667mkdRQW8JQ8vhp4pp/XF4Sg0R/PIdiUl5eTm5tLbm4uDoeDvLw8bDYbaWlpnjd/N9516juSkZHBypUryc3NpbKykpKSEp/HuUsYZ2VlYbPZ2sVw3Ovp6emdnt8ZXZ3TVVs9pasyz256eg86u+fH6+CBPg9FNidLVnp9/0WH3GOq72YJgtAfysrKPN6IxWLxiMeiRYuOiQWUlJSwdOnSY67hDtp7B88BT0DcarXicDjaneOOxXhvd3sW7ge/976uusTc57hHwrlxz1Hpqq2u8GV3Z3R3D7zp7J4fr/THc0kEvId6XAuUcdRb0cecIQhCQKisrKSoqAiHw+GZOW6z2TwPw7S0NM8bdGlpKQUFBZ43//z8fNLT06msrKS0tBSbzeZ5o09PT/fENtxLcXExixcv9ngV7v1uSktL240I87UvJyfH84BfvXo1GRkZncYoKioqKCgoICcnB2jvnXTWlvf9cMeZioqK2LhxIwsWLMBut7ezu7y8nJKSEhwOB2lpaR5bsrOzO70H2dnZ7e5XV/f8eERp3XcNUEpd2FnySqXUNb66zIKBUsoK2LXWnb6uZGZm6o0bN/b62rqlxVgJC4PwcJQ6fh027XSinU7jHkREHNf3QhCOF5RSFVrrzI7b+xtzSeyksZnAoBiHp5TKBgqAQqB7n7mXfLUonyPvv+/doCEybrHx/gwLg/AwVFj4MZ9qWAzh8SMIHxFPWPwIwuPjCRsRT3j8CONzhHub+RkfT9jw4cZ1faBdLlxHGnE11OOqq8NVX4+zvgFXvbleV4+rvt7cXoervgFXfT26pQXd1oZ2OqGtzbOu21qhzdnJvjZoa4OOLyru3+7+/RERqPBwzzbCw1DhHbeFoyIiCE9IIDwxkYikRMITkwhPTCQ8KZGIJHM9MZHwhIROf78gCMGlW3FRSo3obNSXGbi/BqPmy0YgGWNY8g6t9Wa/WtpHtNblSqm8QF3fcs01DJ99JtrlAqcL7XKC0wUuJ7rjZ4djvD91YxPOujpadu3GWVeHq7YW15EjXTeuFGHx8R6xwenEWX9UNI552PsgLDbWEKq4OMLihhMWFY2KiSYsPAIVEQER4aiISEMgOnxvty8ywhSICNAutNOFdrYZguRymp8ucLahvbc5nWB6PLicxr7WVpw1NbTs3Imzurrz+xAWRrjFYoiOJZHwJC8RSkkhOi2N6KlTiehi6KkgCIGhJ57LUnPxibvry/RWqgPVFaaUygWytNbHZIpTSi3G8JSSTJv6n5WuhyRc/u2AXVu3tRlCU1eHs7YOV13tsd+9PlVEOGFxhlCEx8cRNjzu6Hrc0SU8Ls4QlNhYQyQGOa7mZpzV1Tjtdtrs1cZ6tZ226mqcdmO7s7qa5h07jH0OB7hcnvPDk5OJnjqV6GlTiZk2zVifMoWw4cOD+KsEYWjTE3HJVUo9pbX+qKuDtNab/GRTO8xurXQgBx9dbUqpQmCD1nqN+7tSKtf9PZRRERFEJCZCos/ex+OGsOhowkaPJnL06B4dr51O2qqqaNmxg+atW2naupXmbdtxlK5BNzZ6joscN47oadOInjaV6KmG8ERNmoSKjAzUTxGE44aeiEsKUGkGZ8sxRoSV++r26qoLra9orcuBcqVUMuAroc+iDt5MGUaMJeTFRegbKjycyJEjiRw5kuFz5ni2a5eL1r17ad661Vi2baNp61bq33wTnE7joMhIoidPNsRm+skMP/dcoqdOlcEJgtBLeiIu6UA2YAeygJsxUulrjopNmenZFAK3BMjWY1BK+RpIbsewVxDaocLCiBo/nqjx44mfN8+z3dXSQsvOnabobKN561aObKqk9sUX4d77iBw/nvgLLyDugguJzUgXz0YQekC34qK13gmsNAP3T2mtlyilEjC6qbI5KjbVGJ7FgIkLRoylYy4GB4BSyqK1dpixmkxzm11r7TvPhXDcEhYVRcyJJxJz4onttrcePEj9G29Qv+41qlc9hf3RxwhLSCDuvPOIv/AChp97LuFxcUGyWhAGN72a52ImqZxJh8SUXmJTpLUOyNAcM7Zi0Vrne23LBVZqrRO9tlkwRq+laa17PBy6r/NcStZ/wpgTkjjv5LG9PlcIHVxHjtDw7rvUrXuN+jfewFldDZGRDM/KIu7CC4m/8AIiU1ODbaYgDDidzXPp0yRKpdQ8QHecQKmU+rvW+ua+m9llm77EJRso7SAuVozMAYldTZrsSGpqqt6/f7/n+7Jly7qtSulyuTjplr/RbJnE2YkNFP30GuJjpMtkqKOdTho/+oi6deuof+11WnbuBCD65JOJv+AC4i68kJhTpkucRjgu8Ku4mBdMwOgW82Q/VkpNNrvR/E4n4pIOVGitVVfbekKfPBetee3Pi7jljXCaTr6MyCYHyy6ZzA2XntO76wghTbNtJ/Wvv07da6/RuGkTuFxEjBpF3IUXMCInh9g5c0RohCGL38XF68I+vRh/40tczO3VHTyXbKBAa53Tm+v3tVsMwHnYxm9+vYTHnBcSnjSeiQ1b+NfSGxif2rOhs8LQoc1up/7N9dS/to76t99BNzYSNSWN5B/8gBGXX05YdHSwTRQEv9KZuPgjd8YOwKKUKjEzJQ80xWbsxU0O0OsiCvv27UMp1W1XmC/CU6z84cHVvHfLVCbteZFdMdM455513PH75bja2np9PSF0iUhKwnLVdxj3178y7f33SC28BxURyf5f/ZrtF87j8IMP0lZdHWwzBSHgdOu5KKVGAFZzyfJat3J03kkNxgTHMq11p7P5+2Sg0c2VDeRjjA5bjjGgwDvdv3uGvhVw9GWGfn88l3ZoTcljRdzxdgvO5DRi923k77lWzrvihv5fWwhJtNYcef99qh5+mIb1b6FiYki46jsk33gjUZMmBds8QegXfe4WU0q5MNLn24Cd5ucO87MyUDGWgcZv4mLS3NRE/l3383rrNHRbK1mHX+ThO29jxJRZfmtDCD2at22j6pFHqH3+3+i2NuLmXUjyTTcxLD1d4jJCSNIfcdmO0c00pMSkI+7RYj0ZJdYbKrds54cPvkr18Im4vv6IxSMruHXZ/6ISZOjy8UzboUPYn3wSx5OrcNbUEHPaaST/103EZ2cbCUIFIUToj7jco7VeYq7PAyabu+wcO9/laq11SJY29rfn4o3Lpbn7iZf5x+YGnCiSP32SR3JHc/oNd0NMQkDaFEID15EjONauxf7oo7Tu/orIsWNJuvH7JFx9DeFxklhTGPwEcihyEka3mQIWa62n9sfQYBFIcXHztb2BH9z/HNtbEmjZ9yWXfvMP/rz4R4yYeytERAW0bWFwo51O6l9/nap/PkxjZSVh8fEkXruAxO99j8hRo4JtniB0SsCGIns1kI4R0A/J4hkDIS5gBHdXvbOVZc99QgsROCtK+f3Yd7nxjv+07EVCAAAgAElEQVRDnXxZwNsXBj+NmzdT9fAj1JWVQVgYCVdcwQm3/6THWaEFYSAJuLiYjQRshn6gCVTMpTPsDS3c/sh63t7TTMuh3Uz+6K88mZ/FxO//DaLjA96+MPhp2bMH+6OP4Vi9GsLCSP6v/yL5Rz8kLDY22KYJgoc+iYtSapLWelc/G+73NQaCgfJcOlK+5QA/e+IDalsVTRXP8ruRr7PwnlWoiXO6P1k4Lmj5ei8H/3gfdf95mYiRIznhZz8j4corpMSzMCjo6yTKQj+07Y9rDFmyp4/mnV9/i29PT2ZYVi53xS1h3vwfsb/kF9DWEmzzhEFA1LixjPvzn5n45BNEjBrF/qVL2ZU3nyNBeBkShJ7SnefyKlDRn+sD6Vrri/pxjQEhWJ6LN29++Q23PfoutW3htFSuZcXED7lh+RoYeVJQ7RIGD9rlovaFFzj4pz/TduAA8RddxMhf/oKo8eODbZpwnNLXbjG/jJPVWtf44zqBZKBjLp1R19RKwZPv89LWWlqr9nD6Fw/w2P/cQMpFPwfpBhFMXI2NVP3zn1Q99A9oayPx+zeQcvPNhMdLvE4YWAYkoB/KDAbPxZs3vjjAbY+8S52OxLn5Oe4/eQu5d5WCTL4UvGj95hsO3f8XatauJTwxkRNu/wmW3FyZiCkMGIFMXCkEgPNPGs17d17OpVPjiZh5Ff/dcCNXXX05Ne8+GmzThEFE5KhRpC6/m0lrSom2Wjlw513svOoq6t9+J9imCcc54rmYDDbPxZvXP9/PbY+8TT0x8PHzPJC5n8vueByGJXZ/snDcoLWmrqyMg/feR+uePQyfex6jFi8mOi0t2KYJQxjxXEKYC04ewwe/vYqLJsegTv8Oi765guuvuZSGT/8TbNOEQYRSihEXXYT1xRcY+ctf0lhRie2KKznwu99Lmn9hwBFxMelPPZeBIC46gpU35/Dw988gISGBd2b8D6cWPMtr9/0AWhuDbZ4wiAiLiiL5h/9F2quvYJmfR/WqVdgu/RaOtWuRngphoJBuMZPB3C3WkYbmNm5f+Srrvta0Vu/jooOP8+Cf7iNmsqTzF46laetWDiy7k8ZNmxh+1lmMvutOGbos+A0ZLdYNoSQubtZ9sofbHnmbIxFxRG55kUcvH87ZP1oBYeHBNk0YZGiXi+qnnuLQH/+Edjo54fbbSfr+DTKqTOg3EnMZgsybMZ6Nd+dxwWgnbadczoLK6dx0/ZU0H9gabNOEQYYKCyPp+uuxvvgCw+fM4eCKFexacC1NW7YE2zRhiNJrz0UpdSFGnXoLRrp9MGq77MCo77LZrxYOEKHouXiz7uNd/OTh9RyJTibCtp6/fyuJ7BsXg1Q3FDqgtabulVc48Ps/4KyuJvm/biLl1lsJi4kJtmlCCNKvbjFzpn4hRqGwMo6WO3aYh1gw6tdbgSygCijwLiQ22Al1cQFobnPy//72DC98HYlubWZ29cs88qffE5ssEy+FY3HW1PDNvfdSs+ZpIidMYMxv72L47NnBNksIMfpTifIaDFFZ2dM0LqYYzQeqQqUy5WBJ/+IPNm/fyw/+8hyO4RMJO/AZf8y2cNX3FgXbLGGQ0vD+B+xf9htad39FwtVXM2rxLwm3WIJtlhAi9DW32DVApdZ6Zx8bTQDmhYLADAXPxRuXS7Psgcd5bEc0hEdzat27rLqngBGJIVnLTQgwrqYmDj/wIFX/+AfhFgujf/U/xF9yCUq6VYVukNFi3TDUxMWNbfcerr/7cQ4kno5yfM1vzh3BTd+9LthmCYOUpi++YP+vfk3Tp58Sd/75jF72GyLHjAm2WcIgRkaLHadYJ47n/aI7uH3cblxhkdz1yQguuG05Bw7LjG3hWGJOOolJT61i5JICGj74ANtl38b+ryfQTmewTRNCDL95LkqpMwBktNjg5cAeG9f++m/sHDkX1VzHTzKG8/9uuFK6PgSftHz9NQfuvIuGt99m2OmnM/p3vyVm2rRgmyUMMgLmuSilZiqltmOMJluhlNqglJrU3+sK/mf0eCtvPHwfvxn7CbruEH/dEsnsn/2dbXurgm2aMAiJGjeO8SuLSV1RSMvu3ey86moO/Pa3kqdM6BH99lyUUsu11ku72zbYOR48F28cuz7luiX38lnq5YSFhXH9yTH89geXEBEuPaXCsbRVV3P4//5G9VNPERYbS8qtPybp+utRUVHBNk0IMoGMufh6Ipf74boDymBPXOlvLJNO5T9P/pM/jn0bvt7Mk9s0GYv/xaadB4NtmjAIiUhMZPSvf4X1ubUMO/10Dt5TiO3yK6h77TVJhin4xB+ey48wxMRubrICmcBq8/vSUPBijjfPxZvare/wg8W/58Nx1xMRM5xbToum4PuXBdssYRBTv34939xTSIvNRuyc2YxasoSYE08MtllCEAjYUGSllB3YAHSMCmtzW4bWetBPrjiexQWA1kae/eNP+NlnaTD2NNLqN/Ps725hRGJS9+cKxyW6tZXqktUc/utfcdbVYcnN5YT/vp2I5EH/313wI4EUl3la63VKqQRfM/jd+/vVyABw3IuLid22iavvfIxdqdmEHdrG3y5QXPrd2yVHmdApToeDQw88QPWTqwiLiSHllptJvOEGwiQec1zQr5iLUmq1UsqplHpZKTXCe5+XcGQqpa7ueG4oCItwlCTrTN547M/cnPoVzoSxLPrwBH543eW07BbhFXwTbrEw+o47sD7/HLEZGRy89z5sl32b2rIyicccx3QrLkqpX2IkprwXmEInwXpTRKQPZYiw5PZbeP6W2YwIa6Z8Yj4zFv6Ryj8tgLpvgm2aMEiJtloZX/R3xj/0EGEx0ez9ye189f0bJa3/cUpPEleWaK0XeH1/Ffg7RkbkQnOzzVzStdYXB8jWgCLdYr450tLG9+//NxvtUTTt+JAf1RazbOnPiTz7NoiUFO2Cb3RbG47SUg795X9x1tSQcPVVjPzpT4k44YRgmyb4mf5kRf671vpmr+8JQKn51cbR0WEWYJHW+iG/WT2AiLh0jtaav5dvYUX5DlpqDnHC2/fwxGWKU268F06+QuIxQqc4a2s5/MCD2J94grDISFJ+fAtJN96IiowMtmmCn+hPzMXh/cUdtNdaX6S1vtn8TNJah4WqsAhdo5TilpxTePrWc0k5YRTVl6zg7FcmUnj7Apz/vAz2fxRsE4VBSviIEYxaUkDav58n9swzOXjfH9l59dUckRe5IU9PxMWXa1PqY5swxEmfkMhrBTlkTU7GcslPWRGxkLN//xZf/uEceO42iccInRI1aRLjH3yAcQ/8DWdDA7u/dwP77vgfSSUzhOmJuOQrpX7uTkxpIsmojlOS46IpueVcfnx+GvFnXMLeC+8mvWQEf37wYVx/mQlv/QnamoNtpjBIib/wQtJeeIHkhT+i5vnnsV1yKY41a9AuV7BNE/xMT9O/3AJUKqWqlFIlQI5SamLHg5RSF/rVugHkeEv/0h/CwxSLLzmJ4hsySBg7ldE/+F/u2HYy5z/Rxo7SZfBQNhzeFmwzhUFKWGwsI3/+cyY/8zRRU6aw/1e/Zvd3v0fTl1uDbZrgR3oS0L9Ha71EKWUBsr0WK1CNMTT5VaACI9XLgk4vNoiRgH7f2HW4gZv/VcGXB2pp3PgMde+u4r6LY7k5IxL1rRUw8wYJ+Audol0uap5dy8F778VZV0fSjTdywq0/Jmz48GCbJvSQ/owW62zmfQKwgKNiYwG01jrcPyYPLCIufaexxcn/PPsJz2zay/DaXXzxyB1cNCWaf+Q0kjr7Gvj2/TBMarILndNWXc2hP/0JR+kaIsaMMcosz5sXbLOEHhDwMsdKKSuw2lcjoYCIS//QWvPEB1/x239vIUK3sv/FvxC2630evFgxf84kuHolTJwTbDOFQc6RykoOLLuT5m3biLvgAkb/6n+IHDs22GYJXRDwMsdaaxshmGpf8A9KKb43eyLP/+RspqQmkXDp/yMl9y6ueymS6/+1l+oHL4HXl4OzLdimCoOY2PR0Jj/zNCN/+Usa3n+fHd++nMMrV6JbW4NtmtBL/Oa5hDriufiPNqeLh97eyZ/KtkJbCwde+j9G7H2Ph7+lyTn/HLhmJVgmBNtMYZDTum8fB+6+m/rydURNSWPMsmXEZmUF2yyhA33yXPxRrlhKHh9/RISHcfPcNF66/VxOnZBC0rd+StTly7j0uWHctvIdjvzlLPj0mWCbKQxyIlNTGf9//8e4Bx5AH2lk9w3fZ9/SO2g9KAXtQoEuPRel1Exgsta6T08CpdQ1wA6t9eY+2jdgiOcSGJwuzSPv7uLel7+grbWFb17+G2O+eYd/Xa6Y9e0b4dIVECUjg4SucR05wuEHH6Tq4UcAiL/gAiwLFjD8rDmoMCnNHUz6NVoMWApsxwjY13Zz/AiMUWRWoEhrvauvRg8kIi6BZXdVA4vXfMwHO+3o/Vs48PwfWXKag19dNZ3I+f+E1DO6v8jxRFMN7H4Pdr0FNV9D3EiIGwXxYyB+FMSNhvjRMCwJjqOHa8vu3VSXrKbmmWdwOhxEjh+PJS8Py9VXEZGSEmzzjkv6PVpMKTUPyMNIq68xkla6Z+onY4hJIsbclyKt9Wt+sLvfKKVyMfKjWYFyc+DBMYi4BB6XS/PEh1+x/KXPaW5q4mBZMdMOv8G/vhPFydf9Dmbfelw9KNvRVANfvW+Iya63jXxt2gXh0WAZDw2HjGM6EhZpio5bcEwBihtliI9bkIanQFhIzhLwiaulhbpXy3CUlHBkwwaIiCB+3jwSF8wndvZs8WYGEL8ORVZKTcZ4WLsnLzgAm9Z6Z7+s9DPm8Oh8rXWB+b1Ua53n61gRl4Hj6+ojLH3mE97adpi2fVtw/Ocv/CHjMD+57hLCri4yHpBDnabaDmKy2RSTKBiXBZPOhUnnGOvu0gatjVB3AOq/MT7rDkD9ASOnW/2Bo9sa7ce2Fx4NSVZIToOUqZA8FZKnGOuxoV2Gqdm2E8fq1dQ8+yzOmhoiJ0wgcX4eCVddJSWXB4CAz3MZjCilFgMOrXWx+X2H1jrN17EiLgOL1prSjV/z239/RkNTE1WvP0ym/RUeWTCKCTc9BFNzgm2if2muay8m+zaDdhqex7gsQ0gmn2uKybD+tdXWDPUHvcTnADi+gqodULUN7DvB5TW0d1iiITYpUw3xca8nTg6pmj2u5mbqXn0VR8lqI+tyZCTx2fNIXLCA2DPPREmmiIDgF3FRSo3oGHPxtS0QmN1bWW4vpMO+xRjddEkAXmJSCGzQWq8xv+8AMrTWjo7XEHEJDgdqmrjj2U947YuDtO7/giNl/8v/nnmQ711xISrzJjj5coiIDraZvad2H+zbBHs+NASlnZhkGmIyyRSTqNiBtc3ZBo7dULXdWA5vO7pet9/rQGUMGXd7OKNOhYlnGR7QIH9QN+/YgWP1ahxrn8NVU0PUxIlY5s8n4arvEJEU2p7aYKM/Af0N5moJUNkxlmLGYnKAXMCutZ7lH5M9188G0s02bFrr/A77OwqI57u5vsPbc0HEZdChtea5zfv4zdpPqG1spnr945x2+FVWzHUy58RRcMb1kP4DSJkSbFN9c8QO+yph7ybzs9LwGADCImCsW0zOgfFnDryY9IbmOlNodpiis+3o95Z645j4MTDxbJh0tvGZMm3Qio2ruZm6V16humQ1jRUVqMhI4nNySLrpJobNODXY5g0J+iMuLsDSg1Fi6RgP9YBEDU2hsPgQl2qtdaLX92ygQGud46NbrN2x3oi4BJ+DdU38eu2nvPLZN7gO7mD/2kIuSwvj7jPrmJ6ijTf9jB8E15tprjOC7XsrjwqJY/fR/clTYWw6pM6E1HQYPWNwi0lP0RoOfQm73zGWXe8cFdDhJxgezcRzjM+R0wflwIzmbduoLi2lZu1zuGpribvgAk74yW3ETJ8ebNNCmv6Iy2qt9fweNrIxULnFfImLKWjrOohLOlChtVZmQL9Aa51vZnVeKQH9wc/zH+3jV89+QkNTC7XrH8X+wVp+cEkGd2Y6GK8OQGwynPFdQ2iSfYbQ/ENrE3zzaXshObwVT/28hAkw1hSR1JnGcOqYhMDZM5jQGuy2o0Kz+x2o2WPsG5YIE8466tmMnjGoRqo56+up/te/qPrnw7hqa4nPySblttuIOfHEYJsWkvRHXO7RWi/pYSM9FqLe0om4ZGMMe07z2mYFdgCJWmuHUmoRRjwmHVgjQ5FDg29qmyh4+mPe+PIQI11VfPKPJegGO7ffcCVLznSStPc1I4Yx+TzIuAlO+jZERPWtseY6owvo8Fbj7fzwVmOx28Bl5kKLG3VURNyeyXCZV9GO6t2w+13Y/bYhONXm4NHoETBhtiE042cZsZuYEcG1FXDW1WF/7DHsjzyKq66O+IsvJuXWHxMzbVqwTQsp+iMuv9Ra39vDRv6utb65jzZ2d21f4pKL4Y14ey4WjLk2aZ0JiS9EXAYfWmtWfbiH37+4hTA04w6+x6vFvychIYGCn/6Y2+fEEvvZU1DzFcSmwMzvQvqNvr0ZrY0hvB4B2QaHzc/avUePC4swAtYp0+CEkwxvJDUdRqQO2rjCoKV2nyE2u942PJvDXsXAEicbHs3o02DMacZ6/Jig3GNnbS32Rx7F/uijuI4cYcSll5By661EpwXQKx5C9EdcfqG1vs/r+whgJUZwfWmHYx/UWt/iJ5s72tGZ51LaQVzaeS49vb6Iy+Dlq6oj/Lx0Mxt2VTNn/DAcZQ/y8nNrSE1N5c7f/IabzptExEePwZf/Mb2ZuXDqNdBYfdQLObQVmr0mIUbFGQKSMg1OmAYpJxrrSZMhPDJ4P3YoU3/ImM+z/yM48Akc+NjwDt3EphgiM+Y0Q3RGn2a8KAxQl5rT4aDq4UewP/44urGREd/+Nik/voXoyZMHpP1QxS+ei1LqauAhjIf3ImAeUOwO9gdBXDzxla629YTU1FS9f//RYZjLli2TkseDCKdL84+3bdz3ylZGDIvgeyeGs/rPv+K9997jxBNP5A9/+ANXZ89GbX4CKh872v8fN8oUkBOPiknKNPFEBgvNdXDgU1NsTNE5+Dk4W4z9EcNg1ClHvZvRpxvfAzj/pq26Gvs//4n9X0+gm5tJuPxyUn58C1ETj6nsLtBPzwUoxhCVa4Al3t1kSqmFGKPENg90t5i5vdPRYr25vnguocGXB+r4f6s389m+WnIzxpIZtpvf/nopn3/+ObNmzaKwsJDzzzvX6PoakSoVMEMRZ6vx93N7N+5Pd/qb8GhjrtDEs4xl3CyIjvO7GW1VVVQ99A+qV61Ct7aScOWVpNxyM1Hjx/u9rVCmP+LyKpAB7ATyfKV4Mee6WIDsgfRcvLb7nOfSm+uLuIQOLW0u/vraNh54YwejR8Rwz1WnsO2dF1m2bBlff/01l1xyCXfffTczZ84MtqmCv9DayDKw/yPY84ERy9n/kdENqsKN2NjEs4xBAxNmGyPW/ETboUNUPfQQ1aueQrtcWK66ipSb86VCpkl/57ks9o67dHLcZIz4h1+HIpvdXNlAPsYM/OUYCSgrvY5xz9C34jWvpTe4u8WkOyx02PRVNT9f/RG2ww384KxJ3H7+RP5R9CDLly+nurqasWPHMmvWLM4880xmzZpFZmYm8fHxwTZb8BfNdUYGhN3vGsvejWZ3mjK6ztyezYSz/JKvrvWbg1StXImjpAQNDJ89m7i5c4k7fy5R48b1+/qhSr9S7mutfaRjHVqI5xKaNLY4KXz5Cx55dxfWE4bzp/lnMCkeHn/8cd5//30+/PBDtm/fDhilmKdPn+4RmzPPPJNTTz2ViIiIIP8KwS+0NsHeClNs3jGEp7XB2Jc85ahnM/GsflVCbT1wAPtjj1O/bh0tu40JtFFpaYbQzJ1LbPpMVOTxMyjkuExc2RtEXEKbd7Yf5pelH/FNXTO3np/GbRdOJSrCmCVeVVXFhg0b+OCDD/jwww/54IMPqKoyqkUMGzaMjIyMdh7OxIkTJcnhUMDZCvs/NrMKvAtfvXs0bpOaDrMWwilX92twQMuuXdS/+Sb1b75Jw4aN0NpKWHw8w88+2xCbc88Z8nVmAi4uSqlfYgT7QzLHtXSLhT61Ta3c9fwWnq78mlNSR/DgdzOYkHxs6hWtNTt37uSDDz7wCE5lZSXNzc0AjBw5klmzZjFz5kwSExOJi4vrdons4Ztqa2srVVVVXS6HDx9u972lpYXExESSkpJISkpqt95x8d43bFg/sysPNVwuOPQ57HjdGFF4+Euj2Fr6DZD5Q0js32gwZ30DR95/zxSb9bQdPAhKEXPqqR6vJuaU6UOu1kzAxMUcLXYPRsB/ZqByiwUa8VyGDq98doDFaz5GKXjgu+mcldb9m2NLSwuffPJJO+/miy++6HGbUVFRPkUnIiICu93uEY26urpOrxEdHU1ycjLJycmkpKR41qOjo6mursZut2O329utt7W1dXq9mJgYj9BYLBYSEhJISEhot97VEhcXR9gQexB60Bp2rocNK+GLl4xaOtMuhqyFkHZhv3Ojaa1p/vxzQ2jeeJPGjz8GrQk/IYW4c88jbu5chp99FuFx/h/lNtD4XVxMUSnEmPOyRGu9TinlFHERBgO7qxr44aMb2XW4gWVXnMINs3v/VtrW1kZDQwP19fV9XlpaWkhKSvIpGh2X2NjYXnXHaa2pr6/3CE1H4em4vaampt3SlTCBEaPyFpvx48eTmZnpWUaPHt3rezooqdkLFQ9DxaPQcNDI0JD1IyMbt59GnbXZ7TS89ZYhNm+/g6u2FiIjGZGTjeXaa4nNygrZrli/iYuXqGwAVmit13ntC1lxkW6xoUddUyv//dRmXvviIN+bPYFll59CZPgQfRPvJVprGhsbjxEc78XhcLT7vn37drZs2YL7mTFu3Lh2YpORkUFKKMcX2lrg8+fhw5Ww531jAudpeYbQjDndb83otjYaN2+m9tVXPRmao6akkXjtdSRceQXhITaisd/i0kFUlmitN/k4JmTFRTyXoYnTpVnxyhcUvWljjjWZB76bTuLwPia4FKivr2fz5s1s3LjRs3z55Zee/ZMmTTpGcCyWEJzIuv9j2PAQfFIKrUeMiZqzFsH0K/xa7sHV2EjtS/+hetUqmj79FBUbS8Lll5N43bXEnHSS39oJJP0ZinwNhqhUYIjKMZMovY4VcREGJc9Ufs2SZz5h9IgYHroxk2mjQuvtcDBTU1NDZWVlO8Gx2Y7mDJsyZYpHbE455RQmTpzIxIkTiY0NgTo3jQ7Y/KQhNPYdRu2a9Bsh8yZI8O/clsZPPqF61VPUvvgiurmZYTNnknj9dcRffDFhUYP3hahfucUwJi0e46n4OFbERRi0bPqqmkWPV9DY4uQv157BvJP7P7FO8I3dbqeioqKd4Hz11Vftjhk5ciQTJ05k0qRJTJo06Zj1uMEU7Ha5wPa6ITJbXza2jc0wCthNPs+vFUadDgeOZ9fieOopWnbvJjwxEUvuNVgWLBiUkzX71S2mlEoAEjHKGHdakTKUxUViLscH+2saWfRYBZ/uq2HxxSdx81xryAZSQ42DBw+ybds2du/eza5du9i1a5dnfffu3Z6h4G6Sk5N9is+MGTOYNGlS8P5u1bsNb8b2ujFp09UG4VEwLuuo2IzL7Hf3mXa5aHjvPRxPPUXdutdAa4afdy6J115L3HnnocIHx6PWLwF9U2SsQLXWepeP/SErLuK5HD80tjhZ/PTH/PujfVw1cyzLr55BTGRI/rMdMrhcLg4ePHiM6HivNzY2eo4fPXo0c+bM4ayzzuKss84iPT2dmJjAZUrulOZ6+Op92Pkm7HrLzHfmMgYDTDjTEJpJ5xnF5cL7ngmidf9+HKWlVJeW4jx0mMjUVCzXXovlmquJSA7u1EK/DkU2RSYT2OEtMiIuQqigteZvr2/nvle3cvp4C8U3ZDBqRBAeTkKP0Fpz6NAhdu7cSUVFBe+99x7vvvuuJ7YTGRlJRkaGR3DmzJnD2GAklmx0GNkAdq43loOfGduj4mHiHENsJp8Ho2b0aS6Nbm2lbt06qlc9xZEPPoDwcKKtVmKmn0zM9OlEn3wyMSefPKAjzgIyidJMVmkFqsyU+yIuQkjxymcH+FnJZuJjIlj5/UxOGxeCI5uOY7755huP0Lz33nts2LDB0702YcKEdt7N6aef3uNMCn6j4bDh0ex8yxCbqm3G9hgLTDoHTroMTvxWn0pDNG/fTu1LL9H42Wc0b/mctkOHPPsiJ0wgZvp0Yk42RCdm+skB83ACmv7FFJl0jKzIITmRQMTl+OXz/bX86NGNHK5vZkXuaVx5hqRSD1VaWlrYvHmzR2zeffddvv76a8DII5eVlcWcOXO47LLLOPvsswc+A0HtflNs1htpaGq/hrBISLsApn8HTvpWnyduth06RNPnn9O0ZQtNW4zPVvO3A0SMGnVUbE4xhCdizJh+x64GJHGlUmpyV0OVBzMS0D++qapv5pYnKvlwp51bL0jj5zknEhYmgf6hwJ49e9p5N5WVlbS1tZGamkpubi7z589nzpw5Ay80WsPeSvjsGdjyPNR85TehceOsqaHp8y+Ois7nW2ix7TRGvwHhFgsx009m9G9/R9S4vr1USVbkbhDPRWhpc7Hs+U9Z9eEecqaP4s8LziAuWtLxDzXq6+t54YUXWL16NS+99BLNzc2MGzeOvLw85s+fz5lnnjnwI9HcQrPlWfjsuaNCYz0fTvmO0XUWm+SXplyNjTR/+SWNW7bQ/PnnNH22hQmPPdrnPGciLt0g4iKAETh+9N1d/O7Fz5lyQhwP3ZjJ+KQQmOwn9Ina2lpeeOEFSkpKePnll2lpaWHChAnk5eWxYMECMjMzgyM0+yrhs7WwZa1RgTMswhCa6d8x4jR+Ehp/IOLSDSIugjdvbzvMj5+oICI8jKIbMsiaNHj+MwuBoaamhueff57Vq1fzyiuv0NrayqRJk5g/fz7z588nPT09SD/tQRgAABLySURBVEKzCT57tr3QTJ5reDQnfTvoQiPi0g0iLkJHbIfq+eGjG9lb3cjdV88gN2PwzY4WAoPD4eC5556jpKSEsrIy2traSEtL8wjN6aefHjyh2bLW8GocuwEFKdOMxJqpZ8CYM2D0DIgZMWBmibh0g4iL4IuaI638+MkK3tleRf5cK4svPolwCfQfV9jtdtauXcvq1aspLy/H6XQyZcoUzj77bE499VRmzJjBqaeeSmpq6sAJjtawfzNsfcUQnP0fQd3+o/uTpxiCM+YM8/P0Pg137gkiLt0go8WEzmh1urjr35/xr/e/IvvkUdx/rQT6j1cOHz7M2rVrefbZZ9m0aRP79x99oFsslnZi416Skgao26ruG0Nk9m82PvdtNoY6u0mc7OXhmMLjhy41EZduEM9F6I5H393FXf/+jGmj4nnoxkzGJUqg/3inqqqKTz/91LN88sknfPrpp9TU1HiOSU1N9QiNW3imT58+MFmhGw63F5v9m424jZuECZB6Oly8HCzj+9SEiEs3iLgIPWH91kPc+mQl0RFhFN2QScZE/1QqFIYOWmv27t3bTmw+/fRTtmzZQlNTE2BU+ZwyZQqXX3451113HRkZGQPXpXbEbno4Xl7Owtf73G0m4tINIi5CT9l+sJ4fPrqB/Y4m7rlmBlenS6Bf6B6n08mOHTs8YvPhhx/y6quv0traytSpU7nuuuu47rrrOClEioS5EXHpBhEXoTdUN7RwyxMVvG+z8+Pz0/jFRTKjX+g91dXVPP3006xatYrXX38drTUzZ87k+uuvZ8GCBYwf37euqoFExKUbRFyE3tLqdPGb54wZ/ReZM/qHS6Bf6CP79u1j9erVPPnkk2zYsAGA8847j+uuu47c3FxSUlKCbKFvRFy6QcRF6Ataax5+Zxe/f3ELJ40ewUM3ZpJqGRZss4QQZ9u2bTz11FM8+eSTfPHFF0RERHDxxRdz3XXXceWVVw6qKp0iLt0g4iL0h9e/PMjtT24iOjKcld/PYOYECfQL/UdrzUcffcSTTz7JU089xZ49exg2bBhXXHEF119/PZdccglRUVFBtVHEpRtknovQX7Z9U8cPH93Igdom7pXU/YKfcblcvPPOO6xatYrVq1dTVVVFQkICc+fO5bzzzmPu3LmcccYZREQMbNesiEs3iOci+AN7Qws3P17Bh7vs/OTCKfwse5oE+gW/09raSnl5OU8//TRvvvkm27dvByA+Pp5zzjnHIziZmZkBL5Am4tINIi6Cv2hpc/GrtZ+weuPXfGvGaO5fMJOoiJCsoSeECHv37mX9+vWsX7+eN998k88//xyA2NhYzjrrLObOncvcuXOZNWsW0dHRfm1bxKUbRFwEf6K1ZuVbNu5+6QsuPz2V+xecITnJhAHj4MGD7cTm448/BiA6OprZs2d7xGb27Nn9zhQg4tINIi5CIHjwjR0UvvwF35s9gd9deerAZ9IVBIzkm2+99ZZHbDZt2oTL5SIyMpKsrCwef/xxrFZrn67dmbjIoHxBCCC3nJ+G40gLRettJMZG8fOLTgy2SSFFdUMLtsMN7DrcwM7DDeyqaiAuOoIpI+M8S2rCMIlrdUNSUhJXXnklV155JWDUrnnnnXdYv349b7/9NqNGjfJ7m+K5mIjnIgQKrTVLnv6Eko17+PW3p/PDcyYH26RBRUNzGztN8dhpConNFBLHkVbPceFhinGJw6hrasPe0OLZHhsVTtoJcUwdGUeaKThTR8YxISmWiHCJdQUa8VwEIUgopbj76hnUNrXyuxe2kDAsMuQKjzW1OjlQ04RLazRGORHQaA0uDdpc152tY4jswbpmjxfiXg7WNbdrKzUhhkkpw7lsxhgmpwz3LOOTYok0xaKqvpntB+vZfqiebd/Us+NQPe/Zqnhm017PdaLCw5icMpwppuhMNYVncspwYiLDB+zeHa+IuAjCABAeprj/2jOoe2QjBU9/zIiYCC46ZXSwzeqWNqeLpzbs4c9lW6ny8hb6S/LwKCanDGfutBOYlDIca8pwJqUMZ1LycIZFdf/gT46LJjkumjOtye221zW1suNQA9u+qWP7oXq2f1PPp/tqeOnT/bg7acIUJA2PIjE2iqThxpI4PIrkDtvc25Nio3pkU0e01jS1ujjS0saRFieNrU4aW5zmehttTs2Z1mQShgV2qHCwkG4xE+kWEwaChuY2rn/oAz7fX8ujN81iTlpy9ycFAa01r395kLtf+oLtB+uZNTmJ+ZnjiQw3YhtKKRSgFCiU+Wlu9153H2Mel2iKykA/UJtandgONRiCc7CeQ3XNVDe0YG9owX6kheqGFqqPtODq5HE4LDLcFJtIkoZHYxkWidOlOdLS1k40jrQ4aWp1esSkO2Iiw/jWjDFcmzWBrEmJITngQ0aLdYOIizBQVDe0ML/oPfbXNLFq4WxmjEsItknt2LKvlj+8tIV3tlcxOWU4Sy49iYumjwrJB19vcLo0tY2tHrGpajA+7UdasNcfFSF7QwuOxlYiwhSxUREMiwxnWFQ4sVHG57BI93qE8em1PzYqnJjIcGKjImhudfL8R/t4bvM+6pvbsJ4wnGuzxnN1+jhS4vw7FyWQiLh0g6R/EQaSAzVN5P79XY60OFmdP4cpI4OfiPCb2ib++OqXlFZ8TcKwSP573lS+e+ZEmQAaYI60tPHix/sp2bCHjburiQhT5EwfxbWzJnDOlJRBPz9KxKUbxHMRBpqdhxvI+/u7RIaHseaWsxgbpGzKR1raKF5vo+hNG06X5sazJnLbBVNJiB2asYDBzPaDdTz14R6ervz6/7d397FtnHUcwL+/rbRp1ySuu74scZst3bqWZpvqJIWB0AR1YWgVQq0zhrQJgdT0H94EWs0kQEUITZ4mkPiLpAIJ+GNlSSvG0DaRdCBgvCyJu27py1olnWiSplmaNknXNsvLwx/3OLs6tu9i39ln+/uRoiq+586X+vH97nn7Ha5cn0a1bzmaGgJoatiQt/phhcHFAoML5cPJoXE80fofrClfhrb9D2N1DrtDZucUjvQM4Pm/vIuRySk89sBdiDy6BRtX5+DZ7pTW1MwsOk+N4HDX//CPc6MQAR7ZvAZPNG7Azq3r5mfNeQGDiwUGF8qXrvfG8NSv/4t7167EC/s+ifIy91sM/zw3ip+9chqnL05g+0YffvjYVtTX+F1/X1q8C2PX0dZ9AS92D2B44ibuXLkUe4MBfKVxA2rX5L87lcHFAoML5dNfz4xg3++6UV+zCr/9xg7X1mGcuzSJZ189g9fPjCCwajkij27B7gfvKvrB+mIwMzuHv597H4ffvIBjZ0YwO6ew4x4/wvUBfLFufU5uSpJhcLHA4EL59tJbg/juH97Czi3r8Ksng46uLh+9NoVfdJzF4a4LWLH0dnzzs/fia5+6m4sJC9TI5E209wygrXsA50c/QNnHbsMXtq3HnmAg55MAGFwsMLiQF/z+3+/hRy+dxJ5gNZ4PP5RxzqypmVmcHTYWEL49MI6XTwzhxvQsnvzERnwntBn+O/L79EJyhlIKxy9cxdHYAF4+cRHjN6axtnwZvry9GnuDAdy/vtz1c2BwscDgQl7xy2Pn8POOs/j6p+/Gj3d/3LLL6saHszh1cQInh8bROziO3sEJnL00iRm9IrC8bAk+c9+d+P7n78cmD/TRkzumZmbx+ukRHIkN4m/vjmBmTmFbVQX2BAP40kNVWFPuzmQRBhcLDC7kFUop/PTPp/GbN87je7s249s775vfNnlzGqeGJtA7NIGTg+N4Z3Acfe9fm19Z7r9jKeqqK1FXVaH/rcQG/3KOqZSYy9em8KcTQzgaG8Q7g+O4/TbBI5vXYG8wgJ1b1zraHVrSwUVEagGMKaWupirD4EJeMjen8HT72zgSG8BXd2zE5M1pnByawPnRD+bLrKtYhrqqSmyrrsQD1ZWoq67A+ooyBhK6xdlLkzgaG8Qfjw9ieOImysuWYPeDVdgbrEZ9TfYpZ0o2uIhICEAEQFQp1ZmqHIMLec3M7By+9cJxvNo7jMCq5airMgLItupKbKuqwNrysnyfIhWQ2TmFf/WN4mhsEK/1DuPG9CxqVq/Anu0BPPVwTcbjcCUbXABARFoAtDG4UCG6/uEMVixlAnNyzrWpGbzWO4yjsQG8eX4Mb/zgc1hXkdnNiqef5yIiYQCNSqlIkm0HAPQD8AOAUqo1x6dHlFcMLOS0lcuWIFwfQLg+gMvXplzJDJHXHAIiEtLBYz8AX5LtUQD9Sql2HVQ26UBEREQOcCvlUF6Di1KqUyn1HIBYiiLNSql20+8dMAIRAEBEmkXkQJKfkBPnx+zIxDpQ2vj5Z84TYy66heJTSpkDRxDAMaXUqoTXepRSi5rekOmYi4jAC/8/lD+sA6WNn7+1VGMu3kmtuZAfwFjCa1cBQEQWdKGlorvRGgA06eCUrMzBDM8xb9y8o8r22Ivd3255O+XSlcl0mxcV0+e/mH2symW6vdA+f8DbdQDwdsslDOBQQsvFB+AKgE1KqX4H31/V19ejkFoubp5btsde7P52y9spl65MJtu8WgeK6fNfzD5W5TLdXmifP+CdOuDp2WIpJFvwGM8JntiiydZPenp6nhYR84MsLgKAiAw5/F5OqXLx3LI99mL3t1veTrl0ZTLZ5ub/czaK6fNfzD5W5TLdXmifP+CdOlCT7EUvB5cxLJxB5gOAdCvtM6GUOgjgYDbHsJMFgIqDblVfBVALoNPJVjQVDn7n0/PsmItSKoaFrRc/gJSD8vmiZ6e1wBjboSKmLyiNeqZjK4Bovs+Jco/feWueDS5aa8K6ll0wPlBP0bPQePdaGsIA+ky/J50kQsWN33lree0W07O3QjC+sH4R6YPRzRADAKVURK9bCcPoguhLWPeymPdiFgC6RYZ1YjUSLioi4mPXSGHidcE9eQ0uOojEADyXpkzKbXbo5msQRqtnwZ2GnqnWFQ9aIhIVkXCmQYy8z4E6wYfNFzheF9zn5QF9R+jma6eIrEaSFDMwsgCY71o6YGRRjleq5hT7xdItyiTvyrJOXE4o62erpfBke10ga0UfXNJJsahyDEZXHQA2hUuNjTrRDuMiE193xRuMImPnukDWSjq4wCILgN07UlMWAIjIWHzMiAqSVZ3oF5EeU7fKgr56KniW1wV+562VenDxYWH/ebxS+ZF8IecCuh+WzeXiYFknTK1ZtlqKk506wO+8Ba9PRXZbLrMAUGFgnSDWAQeUenDJWRYAKhisE8Q64ICSDi6FlAWAcoN1glgHnFHSwUUriCwAlFOsE8Q6kCVPpNx3kykLwH4Ydx/PwpQFQJeJr8Stxa0DtlSEWCeIdcB9RR9ciIgo99gtRkREjmNwISIixzG4EBGR4xhciIjIcQwuRETkOAYXIg8QEZ+IdIhIm36Ucr7P54A+lwP5PhcqTAwuRDboi39URK7on6j+aXHwIhxTSjXpzMtREenTT2dNd15REVH6X8eCkn5I3z4Am5w6JpWWUs+KTGSLzikV0au2WxKfkKpbHY1KqSaH3i8iIl0AoiIStEjp3p/sMb1E+cSWC5FN+uFgtUieYyoGIJzk9Wy1w1hFnux8QjCekEjkOQwuRPaFYKQBSdaKCMOdxIYtAB534bhErmK3GJF9u5AQQHRrpg1GFl1HusTM9PhLv4iE9QOqzO8bz3tlPp8ggEN6WzzRYi2ATam6zkw5tADAzxxa5AQGFyL7QgA69biLHx8Fkxbzhd8FLTC6xszv0aCU6kwcxFdKxUQkovfpjj9/RERCItKhlNplLi8iUQBd8fPXExfCLv89VALYLUZkg2m8JaqUaldKtZou1K4+QEq3JEL6HOwYgzHIP39eSqlOALV6nAbA/N/UnBBImmG00IiywuBCZM/jMMZb+hNe7wcQzcH7t8O48Me7vrozOEYMQND0ewgfdYfFtQLgzDPKGoMLkT0Lxlu0XC14jHeNAca4iCutJaXUVT7Kl5zA4EJkT6ppvw0w3f3rsY0rItIsImG9uDGYZL9F0d1ayPJYQdwaIBNbMkSOYXAhsqAHzX1I3hUVn7UVX3fSrX869dhMBMZsskwktoraARyKBxoLDeYxGj0JIWaeRq27+FpFpNm8Y8LjfYkywuBClIaeTRUPDs+YB8S1/QCC+gLtM3UpjZnK9C+2xaHfN55eJh4kWmBqeeggEIExUB9NGPDvhhFgQrpc0uwBSqn9AHymlhZnipEj+JhjIoeJSAeAJtM04A4AkXQpXHRgeMaJNC46kEUTpx1ncByfPk7SDAFE6bDlQuQOc5eW3yI3GFHR4SJKInc0iIgfxiwzx1fuE3kdgwuRO17U3WJu5BtLKd4lBmPR5YHE7M1EucLgQuQdQRFpgzE+k7i40Rbd/Zb1Cnudb6wRQFe2x6LSxAF9Igfp2WRtAFr5jBUqZQwuRETkOM4WIyIixzG4EBGR4xhciIjIcQwuRETkOAYXIiJyHIMLERE5jsGFiIgcx+BCRESO+z+JLkgKhn5kKQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x117a99610>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1, 1)\n",
"\n",
"__=ax.loglog()\n",
"\n",
"__=ax.plot(rp_mids, ds_all_sats/1e12, \n",
" label=r'${\\rm all\\ satellites}$', color=mred)\n",
"__=ax.plot(rp_mids, ds_noniso_cens/1e12, \n",
" label=r'${\\rm non-isolated\\ centrals}$', color=morange)\n",
"__=ax.plot(rp_mids, ds_all_cens/1e12, \n",
" label=r'${\\rm all\\ centrals}$', color='k')\n",
"__=ax.plot(rp_mids, ds_iso_cens/1e12, \n",
" label=r'${\\rm isolated\\ centrals}$', color=mblue)\n",
"\n",
"legend = ax.legend()\n",
"\n",
"xlabel = ax.set_xlabel(r'$R_{\\rm p}\\ {\\rm [Mpc]}$')\n",
"ylabel = ax.set_ylabel(r'${\\Delta\\Sigma(R_{\\rm p})}\\ [M_{\\odot}/{\\rm pc}^2]$')\n",
"title = ax.set_title(r'${10 < \\log M_{\\ast} < 10.25}$')\n",
"\n",
"ylim = ax.set_ylim(0.1, 350)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One-halo term lensing for centrals is the same regardless of whether the central is isolated. The lensing one-halo term is a pure mass probe. Not so for the two-halo term. Once you hit a few hundred ${\\rm kpc},$ the isolated centrals begin to show signs of living in an underdense large-scale mode. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [conda root]",
"language": "python",
"name": "conda-root-py"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.14"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment