From b6b998c8c01252fc62929c07e7479c672b479956 Mon Sep 17 00:00:00 2001 From: shouryakhanna Date: Fri, 8 Nov 2024 21:02:14 +0100 Subject: [PATCH] Add files via upload --- .../RedClumpGaiaAllWISE/RC_galfit_demo.ipynb | 90 +++++++++++++++---- 1 file changed, 72 insertions(+), 18 deletions(-) diff --git a/docs/notebooks/RedClumpGaiaAllWISE/RC_galfit_demo.ipynb b/docs/notebooks/RedClumpGaiaAllWISE/RC_galfit_demo.ipynb index cb408fa..036fe69 100644 --- a/docs/notebooks/RedClumpGaiaAllWISE/RC_galfit_demo.ipynb +++ b/docs/notebooks/RedClumpGaiaAllWISE/RC_galfit_demo.ipynb @@ -1,15 +1,19 @@ { "cells": [ { - "cell_type": "code", - "execution_count": 1, - "id": "f8cb6d7d-4f6c-48dc-8ecf-188c39d5372e", + "cell_type": "markdown", + "id": "577d8c67-b5e8-471c-830a-e47de75cec74", + "metadata": {}, + "source": [ + "This notebook demonstrates a use case of modelling the Galactic disc density using Red Clump giants while accounting for the Selection function" + ] + }, + { + "cell_type": "markdown", + "id": "247d2326-b118-4bf9-89a4-16324197f6d9", "metadata": {}, - "outputs": [], "source": [ - "# This notebook demonstrates a use case \n", - "# of modelling the Galactic disc density using Red Clump giants\n", - "# while accounting for the Selection function" + "We begin by importing the necessary dependencies that are placed inside the rcdemo subdirectory. We will also create a temporary directory to place data generated in this notebook, and a directory for all figures generated." ] }, { @@ -22,44 +26,94 @@ "name": "stdout", "output_type": "stream", "text": [ - "1234\n" + "temporary data directory created at:/Users/shouryapro/Documents/pdoc_work/py_scripts/rcdemo/tempdir\n", + "temporary figure directory created at:/Users/shouryapro/Documents/pdoc_work/py_scripts/rcdemo/figdir\n" ] } ], "source": [ - "print('test')\n", - "\n", "import os, sys\n", "loc = os.getcwd()\n", "sys.path.insert(0,loc)\t \n", "from rcdemo.packages_to_import import *\n", - "import rcdemo.dtools as dtools\n" + "import rcdemo.packages_to_import as dtools\n", + "\n", + "# set path to temporary directory inside rcdemo\n", + "tempdir = loc+'/rcdemo/tempdir'\n", + "figdir = loc+'/rcdemo/figdir'\n", + "os.system('rm -rf '+tempdir); os.system('mkdir '+tempdir)\n", + "os.system('rm -rf '+figdir); os.system('mkdir '+figdir)\n", + "print('temporary data directory created at:'+str(tempdir))\n", + "print('temporary figure directory created at:'+str(figdir))" + ] + }, + { + "cell_type": "markdown", + "id": "3e410de2-2322-4f08-b254-fe88f96ea763", + "metadata": {}, + "source": [ + "Let us generate a catalog of mock red clump like stars. We assume here, for simplicity, that the absolute magnitude in RC has a Gaussian like distribution, such that N(Mg,sigma_Mg). Using this, we can generate an artificial Luminosity function which we will later use to sample RC like stars from:\n", + "\n", + "### Luminosity function" ] }, { "cell_type": "code", - "execution_count": 12, - "id": "706ba15d-1e20-44be-a889-e84953000100", + "execution_count": 14, + "id": "6aa0f877-cbed-4399-bff4-e53379e42120", "metadata": {}, "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "lumfunc_use .pkl written to /Users/shouryapro/Documents/pdoc_work/py_scripts/rcdemo/tempdir\n" + ] + }, { "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGwCAYAAACkfh/eAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCE0lEQVR4nO3deXxU9b3/8fdkZ0tCypJAwiIgKCoIBhrsVawgVUP1Pvq7VW/LpnVhsSBUr1jFW22LVUGtIsrPSgqtlfZW8ZJ6UYpVqiJGloooXNkTJIEIJCQkZDu/P/xlSuR8JzPJzDmzvJ6Pxzweyfdz5pzPHIbJZ875Lh7LsiwBAABEkTi3EwAAAAg2ChwAABB1KHAAAEDUocABAABRhwIHAABEHQocAAAQdShwAABA1ElwOwE3NDU16YsvvlCXLl3k8XjcTgcAAPjBsiydPHlSvXr1Ulyc72s0MVngfPHFF8rJyXE7DQAA0AbFxcXKzs72uU1MFjhdunSR9NUJSk1NdTkbAP44VdegUb9YL0n68KdXqmNSTH58ATGtsrJSOTk53r/jvsTkJ0TzbanU1FQKHCBCJNQ1KC65o6Sv/u9S4ACxy5/uJXQyBgAAUYcCBwAARB2u8QKICPFxHn1vRLb3ZwDwhQIHQERITojXou8PczsNABGCW1QAACDqcAUHQESwLEs19Y2SpA6J8UzSCcAnruAAiAg19Y06f8EbOn/BG95CBwBMKHAAAEDUocABAABRhwIHAABEHQocAAAQdShwAESEiy765xw4nTt3Vnx8vBYvXuxiRgDCGQUOgLDn8Xi0d++eFm1NTU2aN2+eOnXq5FJWAMIZ8+AACGvN891YTU2q3vmu9+dmp06dUufOnVVVVeVKfgDCE1dwAIStFpP5Ndar/LVHVP7aI1JjfYvtqquruV0FoAUKHABh6cYbbwxo+3nz5oUoEwCRiAIHQFhatWpVwM8ZP358CDIBEIkocACEHbtCxZOYrL7/Uai+/1EoT2Ky7fP++te/hjo1ABGCAgdA2GlPofLd7343iJkAiFQUOADCypw5c1rdpkePHsbYmjVrgpgNgEhFgQMgrDz77LOtbrN3716fcUZUAaDAARA2SkpKVF9f3/qGkgYMGGCMPfLII8FKCUCEosABEDZmz55tjMXHx7f4/cknnzRue/To0WClBCBCUeAACBtr1641xm6//fYWv+fn5/vcF7epgNhGgQMgbJw6dcoYe/yxx3TF4O66YnB3xf3/GY7/5V/+xbj9Y489FvT8AEQO1qICEBYWLFhgjHXu3FkpifFaPm1Ui/ZFixZp1KhRts8pKysLan4AIgtXcACEhWeeecYYu+WWW2zbc3Nzjc+xLEtFRUXtzgtAZKLAARAWjh8/boz56lDcv39/Y+z+++9vT0oAIhgFDgDXFRYWGmOdOnWSJJ2qa9B5D6zVeQ+s1am6Bm/c162tTZs2BS9JABGFAgeA6xYuXGiMDR482PtzTX2jauobW8SnTp1qfG5lZWW7cwMQmShwALhu586dxtjPfvazVp+fkZFh225Zls+rQwCiFwUOANedOHHCtt3j8bQ6340kDRkyxBj75S9/2da0AEQwChwArioqKlJTU5NtrGvXrn7tY/78+cbY7t2725QXgMhGgQPAVT//+c+NsTP73/ji6ypPVVVVwDkBiHwUOABc9eGHHxpj9913n9/7MfXDqampUUlJScB5AYhsFDgAXFVbW2vbHh8f3+LKTJzHo9H9MzS6f4Z3qYYzZWVlGY9BPxwg9lDgAHCVqYNxly5dWvyekhivVbfnadXteUpJjD9r+5/85CfGY6xZs6ZdOQKIPBQ4AFzja8XvtLS0gPblaz4c+uEAsYcCB4Brnn32WWPsxz/+ccD7M/XDOXnyZMD7AhDZKHAAuMZUeHg8Hs2dO7dF26m6Bo14eJ1GPLyuxVINZzJd9WlsbGTCPyDGUOAAcI1lWbbt3bp1s20/Vl2nY9V1xv2NGjXKGKOjMRBbKHAAuObo0aO27R6bUVL+mDdvnjG2b9++Nu0TQGSiwAHgioKCAmPs6yOo/JWbm6vk5GTbWFuLJgCRiQIHgCsef/xxY2zGjBlt3q9peYfy8vI27xNA5KHAAeCKI0eO2LbbdTAOhOkKTn19PR2NgRhCgQPAFfHxZ0/WJ0mZmZnt2u8111xjjPmadwdAdKHAAeCKL7/80rY9ISHBtj3O49FF2Wm6KDvNdqmGZr7Wr6KjMRA77D9JACCECgsLVV9fbxtLTU21bU9JjNd/z/pWq/vOzs5Wenq67RIQjY2NAeUJIHJxBQeA4xYtWmSMnbnAZlt17tzZtp0lG4DYQYEDwHH79+83xmbNmtXu/Ztucx0/flwlJSXt3j+A8EeBA8BxDQ32Sy107dpV2dnZtrGaukZd+shbuvSRt1RT5/tWU79+/YyxZ555xu88AUQuChwAjjPdKjLdWpIkS5YOnajRoRM1smS/xEMzXzMaM1QciA0UOAAcVVJSYtsBWDIPHQ9Ufn6+EhMTbWOVlZVBOQaA8EaBA8BRK1asMMb69+8ftOP07t3btt3UPwdAdKHAAeCoNWvWGGPtmcH465KSkmzba2pqgnYMAOGLAgeAo0pLS23bk5OTgzJEvJll2ffTKS0tZSQVEAMocAA4yjSCqkePHkE9jukWlcRIKiAWUOAAcNSpU6ds21vrYOyRR4N6dNagHp3lkXmphma+RlK9+eabrT4fQGSjtx0Ax5SUlOjYsWO2MV9DxCWpQ1K81s293O9j5efnKzk5WadPnz4rRj8cIPpxBQeAY55++mljrE+fPkE/XlZWlm17bW1t0I8FILy4XuAsXLhQubm56tKli3r06KHrr79eu3btavV5f/rTnzRkyBClpKTowgsv1Ouvv+5AtgDaw9etoenTpzuWB2tSAdHP9QLnnXfe0cyZM/XBBx9o3bp1qq+v11VXXaXq6mrjc95//33ddNNNuuWWW7R161Zdf/31uv766/XJJ584mDmAQJn636SkpLQ6gqqmrlHjF7+j8YvfaXWphmadOnWybS8vL2ckFRDlXC9w1q5dq6lTp2ro0KEaNmyYCgoKdPDgQW3evNn4nKeeekrf+c53dPfdd+u8887Tww8/rBEjRjAyAghzpqHbOTk5rT9Xlj4/UqXPj1S1ulRDM1+3vVauXOnXPgBEJtcLnK+rqKiQJGVkZBi32bhxo8aNG9eibcKECdq4caPt9qdPn1ZlZWWLBwDneTz2o5+CtUTD182YMcMY8zXhIIDIF1YFTlNTk+bMmaNLL71UF1xwgXG70tJS9ezZs0Vbz549jROILVy4UGlpad6HP98WAQTf7t27bdsbG/275RSo5pFUdsrLy0NyTADhIawKnJkzZ+qTTz7Ryy+/HNT9zp8/XxUVFd5HcXFxUPcPoHUFBQVqamqyjWVnZ4fsuH379rVtD9VVIwDhIWzmwZk1a5YKCwu1YcOGVj/sMjMzVVZW1qKtrKxMmZmZttsnJycbv8UBcMbSpUuNsWCuQfV1pn4/obpqBCA8uH4Fx7IszZo1S6+++qreeustv1YTzsvL0/r161u0rVu3Tnl5eaFKE0A7mUZGdujQIahrUH2dadFNRlEB0c31AmfmzJn63e9+p5deekldunRRaWmpSktLW8w0OnnyZM2fP9/7++zZs7V27VotWrRIO3fu1H/+53/qo48+0qxZs9x4CQDaYeDAgX5t55FHvdM7qHd6B7+WamjWvXt32/aamhoVFhb6vR8AkcX1Amfp0qWqqKjQ2LFjlZWV5X2sWrXKu83Bgwd1+PBh7+9jxozRSy+9pGXLlmnYsGH6r//6L61evdpnx2QA7qqrq7Nt93dW4Q5J8Xrv3m/rvXu/rQ5J/vefGT9+vDHm67YZgMjmsUw3qKNYZWWl0tLSVFFRodTUVLfTAWJCenq6dxqIM11wwQXavn17yI5bUlJiHDk5YsQIn3NuAQgvgfz9dv0KDoDoV1RUZFvcSFKXLl1Ceuzs7GzjmlSmq0oAIh8FDoCQe/75542xoUOH+rWP2vpGffeZd/XdZ95VbX1gI6BMK5XbrTQOIDqEzTBxANFrx44dxthtt93m1z6aLEsfl1R4fw6EaSSVaW0sAJGPKzgAQu7kyZO27enp6crNzQ358U0FzqFDhxguDkQpChwAIWfq62Iawh1spklAJRbdBKIVBQ6AkDNdQXFqhnFfi26uW7fOkRwAOIsCB0DIVVZW2rabCp9gy8/PV4cOHWxjptFdACIbBQ6AkCopKTEucOvkGnFDhgxx7FgA3McoKgAh9fTTTxtj/g4Rb5bRqe1XfEz9gJgLB4hOFDgAQurNN980xvwdIi5JHZMStOUB87ILrTEVMkeOHGnzPgGEL25RAQip+vp62/a0tDRHhog3M032d+TIEYaKA1GIAgdASJk6Eg8YMMDRPEaMGGGMMVQciD4UOAAiQm19o254fqNueH5jwEs1SNLtt99ujDFUHIg+9MEBEFLl5eVB2U+TZWnTvmPenwOVm5urrl276vjx42fFGCoORB+u4AAImXAZIt6sX79+jh8TgDsocACEzIoVK4yxQIeIA0AgKHAAhIyvvi2BDBEHgEBR4AAIGVPfloyMDEeHiDczzYXz2WefOZwJgFCjwAEQMqaConfv3g5n8pVOnTrZttfU1KiwsNDhbACEEgUOgJAxFTi1tbVt2l+HxHh1SIxvcz7Tp083xpYuXdrm/QIIPxQ4AEKmqqrKtr0tI6g6JiXos4e/o88e/o46JrVthoupU6cqMTHRNlZaWtqmfQIITxQ4AEKipKREhw8fto116dLF4Wz+6YILLnDt2ACcQ4EDICQYIg7ATRQ4AEIi2EPEa+sbNW35h5q2/MM2LdUAILawVAOAkAj2EPEmy9Lfdh31/hxs+/fvD/o+AbiHKzgAHOX2cglpaWm27ceOHVNRUZHD2QAIFQocACFhGinlxhpUZxo/frwxtmzZMgczARBKFDgAQsI0UsrNEVSSNHnyZGNsz549DmYCIJQocACExNatW23bT5486XAmLWVnZ2vIkCG2MdMcOQAiDwUOgKArKipSeXm5bSwlJcXhbM6Wnp5u2+528QUgeChwAATd888/b4z56gPjlNOnTwfUDiDyMEwcQNDt2LHDGJs0aVKb9tkxKUH7H7m2rSkBiDFcwQEQdKYrIX369FF2drbD2fhv165dbqcAIEgocAA4plu3bm6nIEnKzMy0ba+urlZhYaHD2QAIBQocABGhtr5RM36/WTN+v7ndSzXMmDHDGFu6dGm79g0gPFDgAAg60wiq9miyLL2+vVSvby9t91IN+fn56tChg22stLS0XfsGEB4ocAAEVUlJiYqLi21jbs9ifCbTXDgAogMFDoCgWrFihTE2dOhQBzMBEMsocAAE1d///ndj7LbbbnMwE9/Cda0sAMFBgQMgqBob7TsAZ2VlKTc31+FszExrYlVUVDicCYBQoMABEFSm5Q769evnbCKtqK+vt23/9NNPVVJS4nA2AIKNAgdAUEXKMggDBgwwxlauXOlgJgBCgQIHQETokBivTx+aoE8fmqAOifHt3t/tt99ujK1bt67d+wfgLgocAEEVijlwJMnj8ahjUoI6JiXI4/G0e3+5ubnq2rWrbYx+OEDko8ABEDSRMgdOs3DrFwQgeChwAATNmjVrjLH2zoFzuqFR8/74D8374z90uqF9SzU0Y6g4EL0ocAAEzYYNG4yx9s6B09hk6c9bSvTnLSVqbGrfUg3NTEPFTe0AIgcFDoCg2b9/v217t27dwmoOnGamIe1btmxxOBMAwUaBAyBoTEPB+/Tp43Am/klJSbFtLy8vV1FRkcPZAAgmChwAMWv8+PHG2LJlyxzMBECwUeAAiFmTJ082xnbs2OFgJgCCjQIHQNCcOHHCtj1cRyVlZ2erf//+bqcBIAQocAAERUlJifbt22cb69Gjh8PZ+G/gwIG27YykAiJbgtsJAIgOK1asMMbOOeecdu+/Q2K8Nt8/zvtzsJhGUpnaAUQGChwAQeFr/aabbrqp3fv3eDz6Rufg3+qKlMVBAQSGW1QAgsK0flNGRkZYzoHTGtOcPgAiAwUOgKAwdSQePHhwUPZ/uqFRD6z+RA+s/iRoSzVIUmZmpm37sWPHmAsHiGAUOACCItTLHjQ2WVr5wQGt/OBA0JZqkKQZM2YYY8yFA0QuChwAQRGpnXXz8/PVqVMn2xhz4QCRiwIHQFDs2rXLtj0SOuuee+65tu2RkDsAexQ4ANqtqKhIx48ft42lpaU5nA0AhEGBs2HDBk2cOFG9evWSx+PR6tWrfW7/9ttvy+PxnPUoLS11JmEAZ/n9739vjPla7ylcmDpIh+sMzABa53qBU11drWHDhmnJkiUBPW/Xrl06fPiw9xHOM6UC0c40g7EkTZo0ycFM2ibUHaQBOM/1if6uvvpqXX311QE/r0ePHkpPT/dr29OnT7e4l15ZWRnw8QCYHTlyxLb9nHPOUXZ2tsPZBM7UEXrLli0OZwIgWFy/gtNWw4cPV1ZWlsaPH6/33nvP57YLFy5UWlqa95GTk+NQlkBsMHXG9fdLiD9SEuL193uu0N/vuUIpCcFbqkGSUlJSbNvLy8uZCweIUBFX4GRlZem5557Tn//8Z/35z39WTk6Oxo4d6/Ob1vz581VRUeF9FBcXO5gxgGCIi/MoJ6OjcjI6Ki7OE9R9++onxFw4QGRy/RZVoAYPHtxiZtQxY8Zoz549euKJJ7Ry5Urb5yQnJ9NZEAihSO+kO3nyZP30pz+1je3Zs8fhbAAEQ8RdwbEzatQo7d692+00gJhlWrcpmJ106xqa9MvXP9MvX/9MdQ1NQduvJGVnZ2vIkCG2scTExKAeC4AzoqLA2bZtm7KystxOA4hJRUVFxmkaEhKCd5G4oalJyzbs1bINe9XQFNwCRzL3Fwr3mZgB2HP9FlVVVVWLqy/79u3Ttm3blJGRoT59+mj+/Pk6dOiQVqxYIUl68skn1b9/fw0dOlS1tbV64YUX9NZbb+nNN9906yUAMe355583xr71rW85mEn7mDpKM5sxEJlcL3A++ugjXXHFFd7f586dK0maMmWKCgoKdPjwYR08eNAbr6ur07x583To0CF17NhRF110kf7617+22AcA5/harykS5sABEJ1cL3DGjh0ryzKvDFxQUNDi93vuuUf33HNPiLMC0F6RMgdOM1OH6DO/YAGIHFHRBweAe0wdiQcOHOhwJu3Tr18/23bmwgEiEwUOgHYxdcKNtM65l112mTHGXDhA5KHAAdAu0dI5d+LEicaYr35GAMKT631wAEQ2Ux+VYE/yl5IQrzfvusz7c7BlZ2erf//+PhcOBRA5uIIDoM2Kior05Zdf2sZMfVraKi7Oo3N7dtG5PbsEfamGZqZ+Q6wqDkQeChwAbeZrDhxffVrCVbT0JwLALSoA7eCrb0p+fn5Qj1XX0KQlf/tqUtCZVwxUUkLwv59FS38iABQ4ANrB9Ie/T58+QZ8Dp6GpSU+t/1ySdPvl5yjJwQvQprW2AIQvblEBaDNTR+LevXs7nElwZGZm2rYfO3aMuXCACEOBA6DNTJ1vI7VT7owZM4wx5sIBIgsFDoA2i7ZOufn5+erUqZNtjLlwgMhCgQOgzXbt2mXbHsmdcs8991zb9kh+TUAsosAB0CZFRUU6fvy4bczUlyUSmPoVBXviQgChRYEDoE18zYEzffp0BzMJrmjrVwTEKoaJA2gTU5+Uzp07B30OHElKTojXazMv9f4cKtHWrwiIVRQ4ANrE1CfF1IelveLjPBqWkx6SfZ/J9Lo+/vjjkB8bQPBwiwoAzmDqP1RdXa3CwkKHswHQVhQ4ANrE6c64dQ1Nev6dPXr+nT2qa2gKyTEk33Ph/Pa3vw3ZcQEEF7eoALSJ051xG5qatPB/dkqSJuX1DdlSDfn5+UpNTVVlZeVZMbs2AOGJKzgA2mTr1q227dHQGff888+3bY+G1wbECgocAAErKipSeXm5bSwlJcXhbIKPVcWByEeBAyBgvubAGT9+vIOZAIA9ChwAAduzZ48xNmnSJAczCQ1mMwYiHwUOgIAlJibatp9//vnKzs52OJvgM3WU/vTTTx3OBEBbUeAACJips21aWprDmYSGqcCpqKhgLhwgQlDgAAiYG51wkxPi9Ydbv6k/3PrNkC7VIEnTpk0zxpYuXRrSYwMIDubBARAR4uM8yhvwDUeOlZ+frw4dOqimpuasWGlpqSM5AGgfruAACNjBgwdt26OpE+6QIUPcTgFAO3AFB0BAioqK9OWXX9rG+vXrF7Lj1jc26Q8fflVY3TSqjxLjQ/v9jJFUQGSjwAEQEF9z4Fx22WUhO259Y5MWvLZDkvR/RmaHvMBxeikKAMHFLSoAAdmxY4cxlp+f72AmoWUaKbZ7926HMwHQFhQ4AAJiGinVp0+fqJgDpzV79+5VSUmJ22kAaAUFDoCg6Natm9spBNXQoUONsZUrVzqYCYC2oMABEJBY6Xx7++23G2Pr1q1zMBMAbUGBAyAgsdL5Njc313hVqra21uFsAASKAgdAQN5//33bdlOn3Eh28cUX27ZHWzEHRKOAh4kfPHhQGzZsUHJysi6++GINHDgwFHkBCEOFhYWqrq62jXXt2jWkx06Kj9OLUy/x/uwEU9EWjcUcEG0CKnB+/etfa+7cuerYsaM8Ho+qqqo0cuRIvfDCC7roootClSOAMPHss88aY9OnTw/psRPi4/TtIT1Deoyvc2PNLQDBEdDXoIcfflj33nuvTpw4oYqKCu3atUvf+ta3lJeXp3fffTdUOQIIE6Z1mDp16hRVc+C0pry83O0UALQioCs4VVVVmjp1quLivqqLBg4cqMWLFysjI0Pz5s3Tpk2bQpIkgPA2ePDgkB+jvrFJq7cekiRdf3HvkM9kLJlHhh08eFAlJSUxMe8PEKkC+oS46KKLtHHjxrPav//97+vjjz8OWlIA8HX1jU26+78+1t3/9bHqG5scOSZz4QCRK6ACZ9GiRZo3b55WrVoly7K87Zs2bdKgQYOCnhyA8BIrc+A08zUXDrflgfAW0C2qb33rWyooKNAdd9yhO++8U8OHD1ddXZ0++eQTvs0AMeDEiRO27dE6bDo3N1eZmZm2fY8aGhpcyAiAvwK+iX3NNdfo888/V0FBgYYPH67ExERJXy2y1717d33729/WnDlzgp0nAJeVlJRo586dtrH6+nqHs3FOv379bNsZKg6Et4DnwZG+uhx9zTXX6JprrvG2FRcXa9u2bdq6dau2bt0atAQBhIcVK1YYYwMGDHAwE2cxVByITG0qcOzk5OQoJydHEydODNYuAYQRX+sv3XbbbQ5mAgCtY6kGAH6pqKiwbc/IyFBubq7D2bhv//79bqcAwIegXcEBEJtMfVSCLSk+Tkv+fYT3Z6ekpaXZth87dkxFRUUxWdwBkYArOAAiQkJ8nK69KEvXXpSlBAcLnPHjxxtjy5YtcywPAIGhwAEAHyZPnmyM7dixw8FMAASCAgeAX0xDxJ3S0Nikv3x8WH/5+LAaHJrJWJKys7OVk5NjG2MkFRC+KHAAtKqwsFA1NTW2sczMTEdyqGts0syXtmjmS1tU52CBI0ndunVz9HgA2o8CB0Crli9fboxNnz7dwUzcEWtLVADRgAIHQKtMs/ampaUpPz/f4WycZ1qKIlqXqACiAQUOgFaZCpzzzz/f4UzcYXr977//vsOZAPAXBQ6AVsX6cgVdu3a1ba+qqlJhYaHD2QDwBwUOALRixowZxtjSpUsdzASAvyhwALSqvLzc7RRclZ+frw4dOtjGSktLHc4GgD9YqgGATyUlJSouLraNOTmKKDE+To/9n4u8PzttyJAh2rp1q+PHBdA2FDgAfFqxYoUxNnToUMfySIyP079dYj/hHgB8neu3qDZs2KCJEyeqV69e8ng8Wr16davPefvttzVixAglJydr4MCBKigoCHmeQKxat26dMXbbbbc5mAkA+M/1Aqe6ulrDhg3TkiVL/Np+3759uvbaa3XFFVdo27ZtmjNnjn70ox/pjTfeCHGmQGyqqKiwbc/IyHB0Je2Gxia9tbNMb+0sc3Sphtbs37/f7RQA2HD9FtXVV1+tq6++2u/tn3vuOfXv31+LFi2SJJ133nl699139cQTT2jChAmhShPA1/Tr18/R49U1Nunmgo8kSZ8+NMHRFcWlryY1tHPs2DEVFRU5WuwBaJ3rV3ACtXHjRo0bN65F24QJE7Rx40bjc06fPq3KysoWDwAIxPjx442xZcuWOZgJAH9EXIFTWlqqnj17tmjr2bOnKisrjYsBLly4UGlpad6HaWVgADCZPHmyMbZjxw4HMwHgj4grcNpi/vz5qqio8D5MQ14BnG3nzp1upxAWsrOzjV+OYmVGZyCSuN4HJ1CZmZkqKytr0VZWVqbU1FTjRFzJycms+gu0QWFhofHKqKlPSjRLTU21ba+rq3M4EwCtibgrOHl5eVq/fn2LtnXr1ikvL8+ljIDo9eyzzxpjvvqkRCtTIcMVHCD8uF7gVFVVadu2bdq2bZukr4aBb9u2TQcPHpT01e2lM+9933HHHdq7d6/uuece7dy5U88++6z++Mc/6q677nIjfSCq+VqGYNKkSQ5mEh6SkpJs248ePepwJgBa4/otqo8++khXXHGF9/e5c+dKkqZMmaKCggIdPnzYW+xIUv/+/fWXv/xFd911l5566illZ2frhRdeYIg44KA+ffooOzvb0WMmxsfpoeuGen92Q/fu3W3bT5w4wVBxIMy4XuCMHTtWlmUZ43azFI8dO5Y1YQAXdevWzfFjJsbHaXJeP8ePe6bx48fr7bffto0tW7aMAgcII67fogKASMFQcSByUOAAMAqnZQgamyxt3POlNu75Uo1N5qu+ocRQcSByUOAAsFVUVKTjx4/bxtwYIn66oVE3/d8PdNP//UCnGxodP34zN27PAQgcBQ4AW88//7wxFotDxAFEFgocALZ89SmJxSHizUxz4ezZs8fhTAD4QoEDwNbJkydt23v37u34EPFwkpiYaNteUVGhoqIih7MBYEKBA8CW6UpFx44dHc4kvFx11VXGGKuKA+GDAgeALdOsvbG+rtudd95pjDFUHAgfFDgAbH3xxRe27abCJ1YwVByIDK7PZAwg/ITbEHFJSoiL0/yrh3h/dlO3bt1UXFzsag4AfKPAAXCWRYsWGWNuDRFPSojT7ZcPcOXYX2fqn2RqB+A8blEBOMv27duNsVgeIt7MVMgcOXLE4UwAmFDgADiLx+Oxbe/Vq5drQ8Qbmyz9o/iE/lF8wrWlGpp17tzZtv3IkSMqKSlxOBsAdihwAJzF1JE4MzPT4Uz+6XRDo65b8p6uW/Keq0s1SNKIESOMsZUrVzqYCQATChwAZzHdgqmvr3c4k/B0++23G2Nr1qxxMBMAJhQ4AM6ye/du2/ba2lqHMwlPubm5xttU5eXlDmcDwA4FDoAWCgsLjfO5xPISDV/Xt29f2/aUlBSHMwFghwIHQAu+hojPnTvXwUwik2W52wEawFcocAC0cOjQIdv25ORk5efnO5xN+DL1Uzpw4IDDmQCwQ4EDoAXTCKpBgwY5nEl4+8Y3vmHbfvLkSVYVB8IABQ4Av5jmxnFKQlycZl85SLOvHOT6Ug2SNHHiRGPsiSeecDATAHbc/5QAEFZMs/G6PYIqKSFOd40/V3eNP1dJCe5/dE2ePNkY+8c//uFgJgDsuP8pASBslJSU6OjRo7axLl26OJxNeMvOzlZWVpZtLD4+3uFsAHwdBQ4Ar6efftoY8zV7rxOamiz9b9lJ/W/ZSTW5vFRDs4yMDNt2RlIB7mM1cQBeb775pjF22223OZjJ2WobGnXVExskSZ8+NEEdk9z/+GLRTSB8cQUHgNepU6ds21NTU5Wbm+twNuHPNOKMRTcB91HgAPAy3Vrp2bOnw5lEhj59+hhjLLoJuIsCB4CXaSg4nWbtzZgxwxhj0U3AXRQ4ALxMi2w2NjY6nElkyM/PV3Jysm2MRTcBd1HgAJAkFRQUqKmpyTbGIptmpkU3ueoFuIsCB4AkaenSpcYYi2yamfotcdULcJf74ywBhIUTJ07YtqekpITFIpsJcXG67bJzvD+HC1O/peLiYoczAXAmChwAksxXInJychzOxF5SQpzuu+Y8t9M4S3p6um17bW2tCgsLw6I4BGJR+HwNAuAqRlC1zfTp042xxYsXO5gJgDNR4ACQJB06dMi2/fTp0w5nYq+pyVLxsVMqPnYqbJZqkKSpU6cai8N9+/Y5nA2AZhQ4AFRUVKTq6mrbWMeOHR3Oxl5tQ6P+5dG/6V8e/ZtqG8KrA69pJJWp8AEQehQ4ALRo0SJjjD4krTMt2WCaIwdA6FHgANDmzZuNsVmzZjmYSWQyLbq5d+9ehzMB0IwCB4AaGhps27/xjW8wyZ8f0tLSbNvr6upUWFjocDYAJAocADLfYunevbvDmUSmOXPmGGOMpALcQYEDQPv377dtD5cRVOGOkVRA+KHAAWJcYWGhsQ9JuIygigS9e/e2bWfJBsAdzGQMxLhnn33WGAunEVTxcR5N+mZf78/hJiHB/uPU1A4gtPifB8S4gwcPGmPhNIIqOSFeD19/gdtpGJn6MR0+fNjhTABI3KICYp5pgr9u3boxgioAptt5zWtSAXAWBQ4Q40xDxDt06OBwJr5ZlqUvq07ry6rTxoVB3XTVVVcZY4ykApxHgQPEONMaVOGmpr5RI3/+V438+V9VUx9+HXfvvPNOY4yRVIDzKHCAGFZQUGC8GtKtWzeHs4ls2dnZSk9Pt40xkgpwHgUOEMOefPJJY+zHP/6xc4lEic6dO9u2V1VVOZwJAAocIIZ9+eWXtu3x8fGaOnWqs8lEAdOQ8OPHj6ukpMThbIDYRoED4Cy9evVyO4WI1K9fP2PsmWeecS4RABQ4QCwzXcFB28ybN88YY6g44CwKHCBGFRUVqaamxjaWkpLicDbRIT8/X4mJibaxyspKh7MBYhszGQMxatGiRcbYJZdc4mAm/omP8+h7I7K9P4ernj170t8GCAMUOECM2rx5szF21113OZiJf5IT4rXo+8PcTqPNysrK3E4BiCncogJiVG1trW17x44dlZub63A20aNLly627XV1dfTDARxEgQPEKFOfkLS0NIcz8Y9lWTpV16BTdQ1huVRDs4kTJxpjv/zlLx3MBIhtFDhADCopKTEWOOFaPNTUN+r8BW/o/AVvhOVSDc1YsgEIDxQ4QAx6+umnjbH+/fs7mEn0yc7OVmpqqm3M4wnfztFAtKHAAWLQmjVrjLH77rvPwUyik6nAYd4hwDkUOEAMOnHihG17YmKi8vPznU0mCpnmEaKjMeCcsClwlixZon79+iklJUWjR4/Whx9+aNy2oKBAHo+nxYOJyQD/NTQ02LabVsNGYEaOHGmMLV682MFMgNgVFgXOqlWrNHfuXD344IPasmWLhg0bpgkTJujIkSPG56Smpurw4cPex4EDBxzMGIhsR48etW2nj0hw+Fqy4fPPP3cwEyB2hUWBs3jxYt16662aNm2azj//fD333HPq2LGjXnzxReNzPB6PMjMzvY+ePXsatz19+rQqKytbPIBYVVBQYIyZ5nBBYHJzc5WUlGQbO336tMPZALHJ9QKnrq5Omzdv1rhx47xtcXFxGjdunDZu3Gh8XlVVlfr27aucnBxdd9112rFjh3HbhQsXKi0tzfvIyckJ6msAIsnjjz9ujM2YMcPBTAIT5/Homgszdc2FmYqLgCtNptt9x48fdzYRIEa5XuCUl5ersbHxrCswPXv2VGlpqe1zBg8erBdffFGvvfaafve736mpqUljxowxrv8yf/58VVRUeB/FxcVBfx1ApDDd+vV4PJo7d67D2fgvJTFez/5gpJ79wUilJMa7nU6rOnToYNve0NBAR2PAAa4XOG2Rl5enyZMna/jw4br88sv1yiuvqHv37nr++edtt09OTlZqamqLBxCrTCuId+/e3eFMots111xjjDGjMRB6rhc43bp1U3x8/FkL0ZWVlSkzM9OvfSQmJuriiy/W7t27Q5EiEDVKSkpUVVVlGwvXGYwjla/5hPisAkLP9QInKSlJI0eO1Pr1671tTU1NWr9+vfLy8vzaR2Njo7Zv366srKxQpQlEBV8zGPfo0cPBTAJ3qq5B/e79i/rd+xedqrMf5h5OsrOz1bFjR9tYY2P4LjUBRIsEtxOQpLlz52rKlCm65JJLNGrUKD355JOqrq7WtGnTJEmTJ09W7969tXDhQknSQw89pG9+85saOHCgTpw4occee0wHDhzQj370IzdfBhD2Vq9ebYz95Cc/cS6RGNG5c2edOnXqrHY6GgOhFxYFzg033KCjR49qwYIFKi0t1fDhw7V27Vpvx+ODBw8qLu6fF5uOHz+uW2+9VaWlperatatGjhyp999/X+eff75bLwGICOXl5bbtHo9HU6dOdTaZGNClSxfbTt2WZWnx4sVh3akbiHQeKwZvvFdWViotLU0VFRV0OEZMMU3kl56eHvZXFU7VNej8BW9Ikj59aII6JoXF9zOfFi9ebJz0r2/fvtq/f7+zCQERLpC/3673wQHgDF9Dk1nqJDR8XaGpqKhwMBMg9lDgADGiuQ+bndGjRzuYSWwxTfhXXV3tbCJAjKHAAWKErzWQfvrTnzqYSWxJTEy0ba+vr1dRUZHD2QCxgwIHiBEnT560bU9MTFRubq7D2QQuzuPRFYO764rB3SNiqYZmgwYNMsZ+8YtfOJgJEFvCv5cegHYrKSlRbW2tbSxSFthMSYzX8mmj3E4jYPPnz9fEiRNtY5s2bXI4GyB2cAUHiAG+lgZggszQys/PN8ZMRSeA9qPAAWLAmjVrjDEm+As9U0fjEydOOJoHEEsocIAYUFlZadseSRP8napr0HkPrNV5D6yNiKUazmQqcCRpwYIFziUCxBAKHCAGmAoc01pJ4aqmvlE19ZG3jtOdd95pjC1btszBTIDYQYEDRDlfVwh8XVlA8Pia8I/5cIDQoMABopyvKwSsheScTp062bZXVVU5nAkQGyhwgCjna40pChzndO3a1RijHw4QfBQ4QBQrKSlRXV2dbcx0RQGhcddddxljS5cudTATIDZQ4ABRzNf8N/S/cRYLbwLOosABotirr75qjEXa7ak4j0ej+2dodP+MiFqq4Uypqam27fX19SopKXE4GyC6eSzLstxOwmmVlZVKS0tTRUWF8QMHiAZJSUmqr6+3jcXgf33X9e/fX/v377eNTZo0SStWrHA2ISDCBPL3mys4QJQqKSkxFjcdOnRwOBtIvufDKSwsdDATIPpR4ABRavbs2cZYRkaGg5mgma/bgr5GuwEIHAUOEKXefPNNYyzS+t9IXy3VMOLhdRrx8LqIW6rhTL5mjy4oKHAuESDKUeAAUcrXBHKRWOBI0rHqOh2rth/2HilGjhxpjC1cuNDBTIDoRoEDRKHFixcbY2lpaQ5mgq9btGiRMbZv3z4HMwGiGwUOEIWeeOIJY2z06NEOZoKvy83NNcYYLg4EDwUOEIV8/ZH8+c9/7mAmsJOZmWmM+eocDsB/FDhAlPE13DguLs7nFQQ44+677zbG1q5d62AmQPSiwAGizH333WeM9e3b18FMYOKrk/epU6cczASIXhQ4QJTZvn27MRbJq1bHeTy6KDtNF2WnRexSDWfyNVw8kv+dgHDBUg0s1YAoUlRUpFGjRhnjMfjfPWx973vf0yuvvGIb69q1q44dO+ZwRkD4Y6kGIEbNmjXLGPPVsRXOe+qpp4wxZjUG2o8CB4giH374oTHmq2MrnJedna2EhARjfM6cOc4lA0QhChwgSrS2WGOkzl7crKauUZc+8pYufeQt1dQ1up1OUIwfP94Y+81vfuNgJkD0ocABooSvKzQ5OTkOZhIaliwdOlGjQydqZCk6+hItW7bMGPO11AaA1lHgAFFi586dxthDDz3kYCbwV3Z2thITE43xG2+80cFsgOhCgQNEgdaGFU+dOtWZRBCwSZMmGWOrVq1yMBMgulDgAFHg17/+tTE2ePBgBzNBoFrra1NQUOBMIkCUocABokBFRYUx9vjjjzuYCdqie/fuxhiT/gFtQ4EDRLjW+mnk5+c7lAna6t577zXGiouLHcwEiB4UOECE89VP44YbbnAwk9DyyKNBPTprUI/O8ijyl2o4U2tD+OlsDASOpRpYqgERbMGCBXr44YeN8Rj87x2xrrjiCr399tvGOP+WAEs1ADHDV3HTp08fBzNBe61cudJnfPHixQ5lAkQHChwgQrX2B2/JkiUOZYJgyM7O9rnC+E9+8hMHswEiHwUOEKHmzZtnjHk8nqjrXFxT16jxi9/R+MXvRM1SDV/X2u1GhowD/qPAASJQa3/o7r//fmcScZAlS58fqdLnR6qiZqmGr2uts/G0adMcygSIfBQ4QARq7Q8dSzNErgceeMBnnL44gH8ocIAI09rEb7Nnz3YoE4RCa8Wpr1uTAP6JAgeIML76aUjSk08+6UwiCJnWitQ5c+Y4kwgQwShwgAhyyy23+Iy3dnsDkaG1IvWpp55yJhEgglHgABHkxRdf9Bmn7030WLRokc94YmKiQ5kAkYkCB4gQHo/v5Qla+4MY6TzyqHd6B/VO7xB1SzXYaW1EVUNDAx2OAR9YqoGlGhABBgwYoL179xrj8fHxamhocDAjOKGoqEijRo3yuU0MfoQjhrFUAxBFioqKfBY3kihuolRubq569erlc5vWruwBsYoCBwhzrX2DZ1h4dDt06FCr2yQlJTmQCRBZKHCAMObPt/NYGRZeW9+o7z7zrr77zLuqrY/OpRpM1qxZ4zNeX1+v8847z6FsgMhAgQOEKX+Km1jqf9FkWfq4pEIfl1SoKYZetyTl5+crIyPD5zY7d+7UJZdc4lBGQPijwAHCkD/FTWvf6hFdvvzyy1a32bx5s77xjW84kA0Q/ihwgDBSWFjoV3GTnZ0ddauFo3X+XLE7duwYHY8BUeAAYaNjx46aOHGiX9sWFxeHOBuEK39vS3o8nlZXnQeiGQUO4LIbb7xRHo9HNTU1fm0fS/1uYM/f98C0adO4moOYRYEDuGT48OHyeDxatWqV38+huEGzQN4LHo9HSUlJKiwsDGFGQHihwAEcNH78eHk8Hnk8Hv3jH/8I6LkUN1JGpyRldGLOl2aBvCfq6+s1ceJEeTwede7cWSUlJSHMDHAfSzWwVANCaNCgQdq9e3e79tGxY0dVV1cHKSNEo/behoqLi9Njjz3W6vpXgNsicqmGJUuWqF+/fkpJSdHo0aP14Ycf+tz+T3/6k4YMGaKUlBRdeOGFev311x3KFPhK8y0mX4/2FjfLly+nuEGr2vs9tampSfPmzfP5XuYWFyKOFQZefvllKykpyXrxxRetHTt2WLfeequVnp5ulZWV2W7/3nvvWfHx8dajjz5qffrpp9b9999vJSYmWtu3b/freBUVFZYkq6KiIpgvw7Isyxo4cKAliQePdj+AQN18882uv2958Gh+eDwe64EHHgjqezyQv99hcYtq9OjRys3N1TPPPCPpq28TOTk5uvPOO3Xvvfeetf0NN9yg6urqFt8mvvnNb2r48OF67rnnWj1eqG5RMVoBwXDzzTfrN7/5jdtphJ3a+kZNefGrK7u/vXmUUhLjXc4ofOXk5NDHBmGjc+fOOnnyZFD2FVG3qOrq6rR582aNGzfO2xYXF6dx48Zp48aNts/ZuHFji+0lacKECcbtT58+rcrKyhaPYBs0aFDQ94nYkpubK8uyKG4MmixLm/Yd06Z9x2JuqYZAFRcXy7Is+hgiLFRVVWnBggWOH9f1Aqe8vFyNjY3q2bNni/aePXuqtLTU9jmlpaUBbb9w4UKlpaV5Hzk5OcFJ/gzt7WuB2HXzzTfLsqxW+50BgaqoqJBlWcrMzHQ7FcS4V1991fFjul7gOGH+/PmqqKjwPkIxC+zAgQODvk9Er3HjxsmyLK7YwBGHDx/2vt/S09PdTgcx6F//9V8dP6brBU63bt0UHx+vsrKyFu1lZWXGbx2ZmZkBbZ+cnKzU1NQWj2D7/PPPg75PRI/Zs2d7/8BYlqV169a5nRJi1PHjx1u8F/v16+d2SohynTt31kMPPeT4cV0vcJKSkjRy5EitX7/e29bU1KT169crLy/P9jl5eXkttpekdevWGbd3imVZXMmJQc23mHw9nnzySbfTBGzt27ev1fcvt7jQFh6PRw888EDQOhgHKsGVo37N3LlzNWXKFF1yySUaNWqUnnzySVVXV2vatGmSpMmTJ6t3795auHChpK++DV9++eVatGiRrr32Wr388sv66KOPtGzZMjdfhiSu5ACIPocPH3Y7BSBgYVHg3HDDDTp69KgWLFig0tJSDR8+XGvXrvV2JD548KDi4v55sWnMmDF66aWXdP/99+u+++7ToEGDtHr1al1wwQVuvQQADujA0HAAfgqLeXCcxlINAABEnoiaBwcAACDYKHAAAEDUocABEBFq6xs1bfmHmrb8Q9XWN7qdDoAwFxadjAGgNU2Wpb/tOur9GQB84QoOAACIOhQ4AAAg6lDgAACAqEOBAwAAog4FDgAAiDoxOYqqefLmyspKlzMB4K9TdQ1qOn1K0lf/dxuSYvLjC4hpzX+3/VmEISaXaigpKVFOTo7baQAAgDYoLi5Wdna2z21issBpamrSF198oS5dusjj8QR135WVlcrJyVFxcTHrXLWCc+U/zpX/OFf+41wFhvPlv1CdK8uydPLkSfXq1avFItx2YvIab1xcXKuVX3ulpqbyH8BPnCv/ca78x7nyH+cqMJwv/4XiXKWlpfm1HZ2MAQBA1KHAAQAAUYcCJ8iSk5P14IMPKjk52e1Uwh7nyn+cK/9xrvzHuQoM58t/4XCuYrKTMQAAiG5cwQEAAFGHAgcAAEQdChwAABB1KHAAAEDUocBpp1/84hcaM2aMOnbsqPT0dL+eY1mWFixYoKysLHXo0EHjxo3T559/HtpEw8SxY8f0gx/8QKmpqUpPT9ctt9yiqqoqn88ZO3asPB5Pi8cdd9zhUMbOWbJkifr166eUlBSNHj1aH374oc/t//SnP2nIkCFKSUnRhRdeqNdff92hTN0XyLkqKCg46/2TkpLiYLbu2bBhgyZOnKhevXrJ4/Fo9erVrT7n7bff1ogRI5ScnKyBAweqoKAg5HmGg0DP1dtvv33W+8rj8ai0tNSZhF20cOFC5ebmqkuXLurRo4euv/567dq1q9XnOf2ZRYHTTnV1dfq3f/s3TZ8+3e/nPProo/r1r3+t5557Tps2bVKnTp00YcIE1dbWhjDT8PCDH/xAO3bs0Lp161RYWKgNGzbotttua/V5t956qw4fPux9PProow5k65xVq1Zp7ty5evDBB7VlyxYNGzZMEyZM0JEjR2y3f//993XTTTfplltu0datW3X99dfr+uuv1yeffOJw5s4L9FxJX82meub758CBAw5m7J7q6moNGzZMS5Ys8Wv7ffv26dprr9UVV1yhbdu2ac6cOfrRj36kN954I8SZui/Qc9Vs165dLd5bPXr0CFGG4eOdd97RzJkz9cEHH2jdunWqr6/XVVddperqauNzXPnMshAUy5cvt9LS0lrdrqmpycrMzLQee+wxb9uJEyes5ORk6w9/+EMIM3Tfp59+akmyioqKvG3/8z//Y3k8HuvQoUPG511++eXW7NmzHcjQPaNGjbJmzpzp/b2xsdHq1auXtXDhQtvtv//971vXXntti7bRo0dbt99+e0jzDAeBnit//29GO0nWq6++6nObe+65xxo6dGiLthtuuMGaMGFCCDMLP/6cq7/97W+WJOv48eOO5BTOjhw5Ykmy3nnnHeM2bnxmcQXHYfv27VNpaanGjRvnbUtLS9Po0aO1ceNGFzMLvY0bNyo9PV2XXHKJt23cuHGKi4vTpk2bfD7397//vbp166YLLrhA8+fP16lTp0KdrmPq6uq0efPmFu+JuLg4jRs3zvie2LhxY4vtJWnChAlR/x5qy7mSpKqqKvXt21c5OTm67rrrtGPHDifSjTix+r5qj+HDhysrK0vjx4/Xe++953Y6rqioqJAkZWRkGLdx470Vk4ttuqn5/mzPnj1btPfs2TPq792Wlpaedfk2ISFBGRkZPl/7v//7v6tv377q1auXPv74Y/3Hf/yHdu3apVdeeSXUKTuivLxcjY2Ntu+JnTt32j6ntLQ0Jt9DbTlXgwcP1osvvqiLLrpIFRUVevzxxzVmzBjt2LEj5IvuRhrT+6qyslI1NTXq0KGDS5mFn6ysLD333HO65JJLdPr0ab3wwgsaO3asNm3apBEjRridnmOampo0Z84cXXrppbrggguM27nxmUWBY+Pee+/Vr371K5/bfPbZZxoyZIhDGYU3f89XW53ZR+fCCy9UVlaWrrzySu3Zs0cDBgxo834RG/Ly8pSXl+f9fcyYMTrvvPP0/PPP6+GHH3YxM0SywYMHa/Dgwd7fx4wZoz179uiJJ57QypUrXczMWTNnztQnn3yid9991+1UzkKBY2PevHmaOnWqz23OOeecNu07MzNTklRWVqasrCxve1lZmYYPH96mfbrN3/OVmZl5VkfQhoYGHTt2zHte/DF69GhJ0u7du6OiwOnWrZvi4+NVVlbWor2srMx4XjIzMwPaPlq05Vx9XWJioi6++GLt3r07FClGNNP7KjU1las3fhg1alRY/qEPlVmzZnkHi7R2NdSNzyz64Njo3r27hgwZ4vORlJTUpn33799fmZmZWr9+vbetsrJSmzZtavEtM5L4e77y8vJ04sQJbd682fvct956S01NTd6ixR/btm2TpBYFYiRLSkrSyJEjW7wnmpqatH79euN7Ii8vr8X2krRu3bqIfQ/5qy3n6usaGxu1ffv2qHn/BFOsvq+CZdu2bTHxvrIsS7NmzdKrr76qt956S/3792/1Oa68t0LWfTlGHDhwwNq6dav1s5/9zOrcubO1detWa+vWrdbJkye92wwePNh65ZVXvL8/8sgjVnp6uvXaa69ZH3/8sXXddddZ/fv3t2pqatx4CY76zne+Y1188cXWpk2brHfffdcaNGiQddNNN3njJSUl1uDBg61NmzZZlmVZu3fvth566CHro48+svbt22e99tpr1jnnnGNddtllbr2EkHj55Zet5ORkq6CgwPr000+t2267zUpPT7dKS0sty7KsSZMmWffee693+/fee89KSEiwHn/8ceuzzz6zHnzwQSsxMdHavn27Wy/BMYGeq5/97GfWG2+8Ye3Zs8favHmzdeONN1opKSnWjh073HoJjjl58qT3M0mStXjxYmvr1q3WgQMHLMuyrHvvvdeaNGmSd/u9e/daHTt2tO6++27rs88+s5YsWWLFx8dba9eudeslOCbQc/XEE09Yq1evtj7//HNr+/bt1uzZs624uDjrr3/9q1svwTHTp0+30tLSrLfffts6fPiw93Hq1CnvNuHwmUWB005TpkyxJJ31+Nvf/ubdRpK1fPly7+9NTU3WAw88YPXs2dNKTk62rrzySmvXrl3OJ++CL7/80rrpppuszp07W6mpqda0adNaFIP79u1rcf4OHjxoXXbZZVZGRoaVnJxsDRw40Lr77rutiooKl15B6Dz99NNWnz59rKSkJGvUqFHWBx984I1dfvnl1pQpU1ps/8c//tE699xzraSkJGvo0KHWX/7yF4czdk8g52rOnDnebXv27Gldc8011pYtW1zI2nnNQ5m//mg+P1OmTLEuv/zys54zfPhwKykpyTrnnHNafHZFs0DP1a9+9StrwIABVkpKipWRkWGNHTvWeuutt9xJ3mF25+nrf+fC4TPL8/+TBQAAiBr0wQEAAFGHAgcAAEQdChwAABB1KHAAAEDUocABAABRhwIHAABEHQocAAAQdShwAABA1KHAAQAAUYcCB0BEmTp1qjwej+64446zYjNnzpTH42l1dXsA0Y8CB0DEycnJ0csvv6yamhpvW21trV566SX16dPHxcwAhAsKHAARZ8SIEcrJydErr7zibXvllVfUp08fXXzxxd62kydP6gc/+IE6deqkrKwsPfHEExo7dqzmzJnjQtYAnESBAyAi3XzzzVq+fLn39xdffFHTpk1rsc3cuXP13nvv6b//+7+1bt06/f3vf9eWLVucThWACyhwAESkH/7wh3r33Xd14MABHThwQO+9955++MMfeuMnT57Ub3/7Wz3++OO68sordcEFF2j58uVqbGx0MWsATklwOwEAaIvu3bvr2muvVUFBgSzL0rXXXqtu3bp543v37lV9fb1GjRrlbUtLS9PgwYPdSBeAwyhwAESsm2++WbNmzZIkLVmyxOVsAIQTblEBiFjf+c53VFdXp/r6ek2YMKFF7JxzzlFiYqKKioq8bRUVFfrf//1fp9ME4AKu4ACIWPHx8frss8+8P5+pS5cumjJliu6++25lZGSoR48eevDBBxUXFyePx+NGugAcxBUcABEtNTVVqamptrHFixcrLy9P+fn5GjdunC699FKdd955SklJcThLAE7zWJZluZ0EADihurpavXv31qJFi3TLLbe4nQ6AEOIWFYCotXXrVu3cuVOjRo1SRUWFHnroIUnSdddd53JmAEKNAgdAVHv88ce1a9cuJSUlaeTIkfr73//eYjg5gOjELSoAABB16GQMAACiDgUOAACIOhQ4AAAg6lDgAACAqEOBAwAAog4FDgAAiDoUOAAAIOpQ4AAAgKjz/wCX5j8anmy2oQAAAABJRU5ErkJggg==", "text/plain": [ - "'/Users/shouryapro/Documents/pdoc_work/py_scripts'" + "
" ] }, - "execution_count": 12, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "loc" + "# Generate RC Luminosity Function\n", + "mgmin,mgmax = -1., 2.\n", + "mg = numpy.linspace(mgmin,mgmax, 10000)\n", + "ye = norm.pdf(mg,loc=0.44,scale=0.17)\t\n", + "\n", + "zm = interp1d(mg,ye,bounds_error=False,fill_value=(ye[0],ye[-1]))\t\t\t\n", + "gridM_G = mg\n", + "gridLF = zm(mg)\n", + "gridLF = gridLF#/gridLF.max()\n", + "\n", + "lumfunc = {}\n", + "lumfunc['gridM_G'] = gridM_G\n", + "lumfunc['gridLF'] = gridLF\n", + "\n", + "# store the luminosity function as a pickle file in tempdir\n", + "dtools.picklewrite(lumfunc,'lumfunc_use',tempdir)\n", + "\n", + "cdfval,mgmin,mgmax = dtools.pdf2cdf(ye,xv=mg,getlims=True,usesig=2)\n", + "\n", + "plt.close('all')\n", + "plt.plot(mg,ye,'k.')\n", + "plt.axvline(0.44,linestyle='--')\n", + "plt.xlabel('Mg')\n", + "plt.ylabel(r'$\\rho$')\n", + "plt.savefig(figdir+'/lfrc.png')\n", + "\n", + "#--------------------------------------------------------------" ] }, { "cell_type": "code", - "execution_count": 30, + "execution_count": null, "id": "1da57cdf-1844-4006-81fd-11b1616b66b1", "metadata": {}, "outputs": [],