Skip to content

Instantly share code, notes, and snippets.

@smsharma
Last active February 27, 2022 19:53
Show Gist options
  • Save smsharma/85ff5df9946415ae2eac8533bee5e2ae to your computer and use it in GitHub Desktop.
Save smsharma/85ff5df9946415ae2eac8533bee5e2ae to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import healpy as hp"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def mod(dividends, divisor):\n",
" \"\"\" Return dividends (array) mod divisor (double)\n",
" \"\"\"\n",
"\n",
" output = np.zeros(len(dividends))\n",
"\n",
" for i in range(len(dividends)): \n",
" output[i] = dividends[i]\n",
" done=False\n",
" while (not done):\n",
" if output[i] >= divisor:\n",
" output[i] -= divisor\n",
" elif output[i] < 0.:\n",
" output[i] += divisor\n",
" else:\n",
" done=True\n",
"\n",
" return output\n",
"\n",
"def make_ring_mask(inner, outer, ring_b, ring_l, nside):\n",
" \"\"\" Masks outside inner < r < outer, of a ring centred at (ring_b,ring_l)\n",
" \"\"\"\n",
" mask_none = np.arange(hp.nside2npix(nside))\n",
" return np.logical_not(\n",
" (np.cos(np.radians(inner)) >=\n",
" np.dot(hp.ang2vec(np.radians(90-ring_b),\n",
" np.radians(ring_l)), hp.pix2vec(nside, mask_none))) *\n",
" (np.dot(hp.ang2vec(np.radians(90-ring_b),\n",
" np.radians(ring_l)), hp.pix2vec(nside, mask_none)) >=\n",
" np.cos(np.radians(outer))))\n",
"\n",
"def rho_NFW(r, gamma=1., r_s=17.):\n",
" \"\"\" Generalized NFW profile\n",
" \"\"\"\n",
" return (r / r_s) ** -gamma * (1 + (r / r_s)) ** (-3 + gamma) \n",
"\n",
"def rGC(s_ary, b_ary, l_ary, rsun=8.):\n",
" \"\"\" Distance to GC as a function of LOS distance, latitude, longitude\n",
" \"\"\"\n",
" return np.sqrt(s_ary ** 2 - 2. * rsun * np.transpose(np.multiply.outer(s_ary, np.cos(b_ary) * np.cos(l_ary))) + rsun ** 2)\n",
"\n",
"def get_NFW2_template(gamma=1.2, r_s=17, rsun=8, nside=128):\n",
" \"\"\" Generate squared-NFW template by doing LOS integral\n",
" :param gamma: NFW spectral index\n",
" :param r_s: MW scale radius in kpc\n",
" :param rsun: Distance from Sun to GC in kpc\n",
" :param nside: HEALPix nside parameter\n",
" \"\"\"\n",
" \n",
" # Restrict to a region around the GC to make computation faster\n",
" mask = make_ring_mask(inner=0, outer=65, ring_b=0, ring_l=0, nside=nside)\n",
" \n",
" mask_restrict = np.where(mask == 0)[0]\n",
"\n",
" # Get lon/lat array within restricted region\n",
" theta_ary, phi_ary = hp.pix2ang(nside, mask_restrict)\n",
"\n",
" b_ary = np.pi / 2. - theta_ary\n",
" l_ary = mod(phi_ary + np.pi, 2. * np.pi) - np.pi\n",
" \n",
" # LOS distance array\n",
" s_ary = np.linspace(0, 30, 500)\n",
"\n",
" # LOS integral of density^2\n",
" int_rho2_temp = np.trapz(rho_NFW(rGC(s_ary, b_ary, l_ary, rsun=rsun), gamma=gamma, r_s=r_s) ** 2, s_ary, axis=1)\n",
"\n",
" int_rho2 = np.zeros(hp.nside2npix(nside))\n",
" int_rho2[~mask] = int_rho2_temp\n",
" \n",
" # Normalize to mean in the region\n",
" int_rho2 /= np.mean(int_rho2)\n",
"\n",
" return int_rho2\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"template = get_NFW2_template(gamma=1.2, r_s=17, rsun=8,)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAFyCAYAAADGe88vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABJKklEQVR4nO29e6w0yXne97zVcy7fbXe5XO6SCnWByUiknBhURNsJHCmWE1KGZScCYoGRHDgInMAKoAg2EBtwYBgIEOQCJxGQwAmc2H8Elh3ZMgTJF1mSBYm6BXEgmbJimCIlWrJEidzlcrn7fd9+33fOTHflj6rqrq6p7ume6Ut1z/MDDnqmL9XVPedMP+d533pLtNYghBBCCCHDoubuACGEEELIGqHIIoQQQggZAYosQgghhJARoMgihBBCCBkBiixCCCGEkBGgyCKEjIKIbEREvPcXI53n/SLy50Xkw2O0Twghx0KRRQipISJ/QET+noh8SkRUsO05EfktEfneNlEjIh8F8GsAvkJE3iUifx/Af3liv75BRH44sumfAxAAv/OU9r3z/GUR+WNDtEUIOW8osgghNbTWPwHg+wB8AsA3BZs/BOBXAPwVrfXPt7TxYwB+wr7+AoD/YYCu/RKA/zZyrgLAbwzQvuMvA/jJAdsjhJwpFFmEkBgawN8A8B3BemW3daliXDS8bkRElB9irHVI67e01j/TpZ1jcK6d1vofa61/+9R2CCGEXwaEkCZ+BMA3icgVAIjIuwF8PtxJRP6kiPxREfkPReTb2hoUkUxE/hcR+TUR+bCI/FUR+SkReUlE/isAPw7gnSLyu0Tku0Xk20Tkz9ljv1tE/pHX1sdE5E+IyB8G8G94698hIn9ORP49EfmvReRBpB8fE5FHIvId9v1fAvC3ReQrReT7ReTPNLVlr/OJiHyLiPxnNnz6IRH5qIh8GsbtI4QQiixCSByt9S2M6PlDdtUHtda/7O8jIt8K4IHW+m9rrb8XwO8Xka9raTMH8KcBPAbwCwD+RwBvaK1fhxF1fxzAmwD+DwD/q9b6+wHcE5EPa63/ZwBbe94PA/io1vqvaq3/HoBf9U7zPwH4Ma31DwH4OIA/GenH34Rx6l63qz4O4GNa638BEy6829SWvc4fB/ArWuv/zfb3/wPwUwD+C631P266fkLIeUGRRQhpww8ZxsJ4HwPwT7z3vwjgj7Y1qLXe2f1+N4CvBiAichfAe7TWnwXwNQBeAPBvi8g3A/gcgMwevrXLfx/AP/KafeS9/sMA3muPfQGVkAr56/61aa1d21tvn6a2fgjAt4rIhwD8LIBvsD8fb7t2Qsh5sZm7A4SQpPlpAH9FRL4edbfIEX6HKABdSjX8EIBvhRFKPwLgm1GJuBzAY631j9r3P7p3tBFdTXlhOYCf0Fo/atju+BkAf0lEPgDg13u29XcB/ACATwH4HgDfCeBntNYPD5yTEHJG0MkihMS4EhFlR+79IIDv1Fq7EXyCShD9AIwj5fh6u/8hfgTGjXodwN8B8GcBfNJu+zSMu/VOmBeviMjvCY7/IQC/y3v/MirB9wPwRkWKyB+JdUBrrWHE0l/QWv+/Df2MtqW1fg1G5D1vQ6gfRLfBAISQM4JOFiGkhoh8BMB3AfjdIvIXAPw1WCElIt8JIyj+YxF5qLX+v0TkT4vIHwdwDeDntNb/t4h8I4B/HcCfEpH/DsCfAPB+EflBrfU/0Vo/FpEfB/D/aK1zEfklrfU/A0xJBhH5GIA/axPdN1rrvyUi/wmAf01E/ojW+u+KyAdsPasvAfhymPDiT8EItj9vE95zxJ0wx1+DcdTctX8ZgP8UwFeJyP95oK2/CZNXBpiQYZNQI4ScKWL+mSOEkMM4d8uVWdD8AiGEkEYosgghhBBCRoA5WYQQQgghI0CRRQghhBAyAhRZhBBCCCEjQJFFCCGEEDICXUs4MDueEEIIIWSf6KT2AOtkEUICPqJa53gmHfiHxffP3QVCSAJ0LeFAJ4uQhUPxlB4UY4SsgkYniyKLkJVAEbU+KMIIWQQUWYQsDYom0heKMkJmgSKLkFShmCJjQ/FFyKhQZBGSAhRUJBUovAgZDIosQsaE4omsFYoxQg5CkUXIkFBUkXOFoouQPSiyCDkFiipC4lB0EUKRRUgrFFGEjANFGDkDKLII8aGoImQeKLrICqHIIucJxRQhy4DiiywYiixyHlBUEbIOKLrIgqDIIuuFwoqQdUPBRRKHIossH4opQogPxRdJBIosskworAghXaDgIjNCkUWWA4UVIeQUKLjIxFBkkTShoCKETAGFFxkRiiySDhRWhJA5oeAiA0ORReaH4ooQkhIUW2QgKLLItFBQEUKWCIUXOQKKLDI+FFaEkDVBwUU6QpFFhoeiihByTlB0kQYosshwUFwRQs4Zii0SQJFFToPCihBC9qHgIqDIIn2hqCKEkP5QdJ0lFFmkGxRXJEnEfodpbV67pVtHSGJQbJ0VFFmkHYorMjvS+D11PBRgZGYots4CiixSh6KKTM4YIupUKMLIxFB0rRKKLGKguCKT01dciTr9nLrouT+/4si0UGytCoqsc4WiikzKIUE1hIAaikNCjMKLTAhF16KhyDpHKLDIZCxJXIVQbJFEoNBaLBRZ5wKFFZmUNnGVsrBqok1wUWyRCaHgWhQUWWuGwopMxtpE1SEoukgCUHAlD0XWGqG4IpMSE1hrFFZNxAQXhRaZEIqtZGkUWWf0DbkuKLDIJDhhlWL5hRTg/SETwu/95UEnayHwj4tMSigaRBkn55ycq0O4+xE6XHS3yITQ3UoCOllLhgKLTEqTK0OBVafpftDVIhPC50Pa0MlKGP7xkMmJOVikG3S0yMzQ1ZoNJr4vBQorMil0rYanaUQiRReZEAquSaHISh2KKzI55z5acGw4GpEkAMXWJFBkpQiFFZkFiqtpodgiiUDBNRpMfE8NCiySDBRY48L7SxKBz53poZM1MfwlJ7PBpPb5YXI8SQS6WoPCcOHcUFyRWWBoME0YQiSJQLE1CAwXzgkFFpkFCqx0iX0OrK9FZoDPp3GhkzUi/OUls+I/tBcmrkQdLzh0sbCvK9/VoptFZoSu1tEwXDgVFFZkVhbkXp0ipI4lWQHG8CFJDAquXlBkTQEFFpmdhJPb5xBVh0hKdDEpniQGhVZnKLLGhgKLzE6CAitFYdVEEoKLQoskBoVWJyiyxoDCiiRBQuJqSaLqELOKLootkiAUXI1wdOHQUGCRJKDAGo1Zryf8HDnykCQAn3v9oZPVE/6SkWRIYPTg2oRVG7M5Wxx9SBKErlYNhgtPheKKJMXMAuucxFXILGKLQoskCsUWAIYLT4MCiyQFBdaszHL9/ufM0CFJCD4f26GT1QJ/eUhSzCiuzl1YtTG5s0VXiyTKGbtadLL6QoFFkmImgSVK0hBY7ppFVT/++hmZ/B7R1SKJwufmPvN/QyUIf1FIUswosCYnFFPhuqZ9m46ZEAotQvj8DGG40IO/HCQ5ZhBYk4iFOR2o2BQ2g59ioq9Mhg5JwpxR+JCjC9uguCJJsjaBlUBob48RBReFFiGGMxBbzMkiZLEsXWD5YbzUGLFvk4UPU723hJDzdbLoXpFkWbqDJcq4K0t8+Lt+D+hw0dEipGKlrhbDhT4UWCRZJhZYo7gtSxRXISOEEScRWxRaZAGsUGhRZAEUVyRxJhRYgztXa2dJzhaFFlkIKxJbzMmiwCJJM+Ew/EULLCXVz5QMeJ3TlnpgeQeSLufwXD4LkXUOHyRZMEt9EI4lsHwR5YuqUJzEto0pwJbq2C3194ucBWt/Pq86XLj2D4+shInChIM5KEP1cY5ip0OF6gYKH44aOmTYkCyMBYcPzy8niwKLLIKlCKwlC6smEhFcFFqEVCxUaDEni5DkmHHC58mZI4/qECn2aWg4/Q4hs7IqJ4vuFVkM4QNvBJE1SHjwlH4tUcCc4ioNEEIcxdUK+0VHiyyEBblajV92myl7QQiJkKKLdW7iyuH6fozYcfdsgrkRezFwcVVCSHdW42TRxSKLYYIw4Uku1rF9GlhcyRHhrY7fZ9051lk6QdSMlqPF/CyyQBbiZq038Z3iiiyKkQXWLCHCI895jIg6laNFWF/hs4TQIYUWWRCJiy0mvhNCWnDCagKBJSKzCKyTzt33Oo+9n4SQVbFYJ4sOFlkUawoR9jzPXIKqD70drj4uE0OHhAxGoo7WupwsCixC6sySg9Wl6VMFllL1ZdNrddo1jCoET7i/k07BQ8gCWNrzf3FO1tJuMCFJl2sYycE6WrQoBRTFaaLJHV8c5yD1crUmcLRY1oGQfRJztNaR+E6BRRZHyonuIwisXuLqRPfpKHoIr85ia8mhQ4YNyYJJSGgtX2RRYJFFMqKLNYnAWou4CplTbKUitOhmkRWQiNBadk4WBRZZJBOECY9iYIHVixQEFjBOP7rer1R/DxYwQIGQkNT1QfJOVuo3kJBG3EMrpTDh1A5WXzEz5IO+94jBdocpFUdr1LAh3SyyUGZ2tJbpZFFgkcWSoiswkMBywupkgVWKUKl+hiRs98T+dr/uhTlaPin+3hLSgVT1QoJ/5YZUbxghc3OUizXwA71VaCjVLlj6CJ+hCAVdEwf6PniphyM+F5Z1ICROirohyXBhijeKkM6MOKJwNIHVwb3SWjeLjC5hwa4CxbXVt5SDv3/XxPYu338Nbbn7cfA7tEt474jQ4eBhQ440JCthhtDhcsKFFFhkNaQQDhpyNOOYAutQ0dEuHFOktIvwa2hjUFdrZb8rhMxJSjoiKScrpRtDyFGMlOw+p4PVfnzDOQZOiG/rR/8pcQ64Rk3tnZoYvzRHi24WWTgTOlrpO1kUWITEmUtgHc0ACebhz5D7D+JuHUOX+80cLUIGIwVdsZm7A4SshpWNzDoqPHgoqbzhPK35Xkfg2mrMmzqUuyWy7+T4/Y8c1yk/a0nE7gEhpBezhwtTUJqEnMxIye6jjiRsabu3wDpCkDWfw64vtHntlk3E9msIobV+3zUJrp7hw/ZzjFdDa7Rq8BRaZAWMHDpMP1xIyGKhwIrv2+JcdXKtXB8P3Yeu+x06d1Ofe96P9vIW49XQGjRs6J9/ZQ4tIVMyq8iii0XIDAz6MG4QWHu7tRTyVFL/GYIDbbb3p4fQOrZvhJBJmUtvzBYupMAiqyClCaBPTHbvLDiqAw7u2yiqGtsc8P++tpBbJLQW/S4MQ4I9w40nhw57hg05gTQh7YwUNkwrXEiBRVaD/+A5F4EVq5retSRDX4HV9T7E9mu7J53bjUyi3CN8eHLosOfv1GghQwosshKm1h+TO1kUWGQ1pJKLNbXAOrBvZ/fqFFHVhZirE3OGurhasST3Hgnxq3G0KLbIShjY0UrLySKEpMGQAuvgaEHACCv3E24PBVZsAulD62Ln7Xr+gNr1nJgQP/ich4SQRTCpyKKLRVZJ6tORdBERPjHXpmeCe3neUOAc6uMh4XRoXXh81/Bk0NeTEuIbSjt0EqEpkvrvNyFHMJUemSxcSIFFVkUKocJD5z2mTAMQz0Nq2x62t+dItbhGTf1ouja/HlZTSKwptLYXAvTeh8d4245KiG+Zguek0GGPsCFDhoQcZqCwIcOFhJCKxgf9IYEVobPACtsMhWrXUN6hdU1ttbpjQT8PuVp7x3cbCLCqivCEkINM4mTRxSKrYwQna/CE91OT3fvmYHV1sA7duy5V3ptoc7ma3JiOjpY5zD+uQzL80EnwcyXA08kiK2YAN6vxy2p0kUWBRVbHOQqsthGEbXlXqsE9CgWYJ6pOmctwb67CmOhqqv8UChJ/v7bw4RGhQwotQtLiRKHFcCEhgzOXwDrlXMc6WE1tdBFYtYObRxb6U96IH9oLw3yxZXBc2VbbCMOQjk7cQfF3bJHWkRjs94sJ8IT0ZlQniy4WWR1rS3jvUq7hGBcr5mC57UH/9vrUNSm/Jbk8dFr23Sf7PubONLleXd2sjvWzmABPSFqc4GZNHy6kwCKrpHRPliWwek36XJY/6Ciu/H7Fcq+CbdFwYMM5IWIe5m3Cy21vEj66Lo72wolAJR5iuVp9Rh5Gzrm3ba/7TaIscaFFkUVWyJFCa9pwIQUWIQvg2ImQDwmsWoivxb3aq2ml6n3oGjJs2ua3Z38ak/Ul2De2ren6Y/Qp8koISYah9ctm0NYIIb2YKkx4uO2OeVhdBFZkfatrFdsWPZ+0LwtdP15rIMvM+9K5KmpHAoB2t7WwL3RRHe+S5kXVnSO7vuaKhbi+HOv4tNUCA/b71IIoGdbNIoR0gv9aEdKVgUOFgwqstkOPcVW6TJdzSGDFO1O1HzpZmdee+3Hb3f5NS5H6cX57ErTj96Ptug45WuUlNeev7TF0EnyP34vBE+AnTN4nZKkMLrIYKiSrZIQHSi9n4ZBj0ecBOqbAivSpNlKw6fx7gkb222wKC/rbY+3s9asuEhpdumi/xhFazfse+FznystyUGiRFTKkjhk08Z0Ci6yWuWtjHREq7JTs3iJCDo4i7BoiDN2r2rqWXCsA2rVVaGgljUu3jzmoYRnbHiSqH0yIjyXDxxLhDxUq7ZMEn2ICPEcZkjOgRxJ845c5c7IISZ0jRV2ngp6HRu419eVAiLAmsJoS2qOJ8nVxVb7P2pfmGHtITHPU8qukum7Xt6IAwpGH7rq65GgdooMQOXrKnb59IYRMBnOyCOlDigVI+5RtOGLy571QWVvflcRDZzXR5eVbRdwtnal94dU0ijCyn85i51T75/b3CUOkB67x4P04dZLt2PkGgIVJCZmWwcKFDBWS1ZJyqPCUulgNeUQHyxzUHKiWEGEooELnKnStlLJukao7bE4Uad289EOBIsaZcksEoUT3A5jjwvBhW+iwKWzYdY7DY0OGkXPUD2LIkJAx6BgybPxC578jhJwLx9RpOpDw7RMNEe7vVBdYoevkjstU9ePWtS39/f1tsfPEkurL662EYadRf13cvSZYN4uQ1TPIXzldLEK6M3bphsHmxYuGweIuVvSctdIJLXlXSlVCSUlNWGlvqUWgNwo6s0uR2vayrbCdcp9IuNEvExG5hlroUDxnru3enMhRn98cpRwIOQNO1TcnhwspsMjqaUjuPqqpvg+4oUKFHXOEoiMKRQWjC+vr90YRugKgERGjVSBoQocJgBap/v0L246xN4oQEB2EBwEgL2rrpND18KEfOszzWpu63N4w/U4kbNg40rBt3sPw2JCBQoZm9wFCfJxih5wJB8KGDBcSchRz1gEay8XqK7Bix3rrdS1/yhNNMYFVayNwtjIxAmtTLzSqlSq37TlZmXWz/Bwwb1vr+UIXy3fbvOvYG20Yu4eRsKGEbYX3sIWx3azBmfPvhJCEYQkHQppYy4OjLfenj8ByhGHCMMk9kv+kfQETbNeZ1F8rew6toTMF0UaASa5rJRtEdLm+EEDyojwOkkG0hr4wifGSa5iG6q6XVqiXfFAAIMb1yjJzrN1fK7SXb3DrvelwypIQSjUmurduWxJuEAIhpOSkf30YKiSrJuZgrAkVuFGN+0VCdrF1/si7MM8p5hBZh8p3sqBQOVg2t6qwzlRxUc+tKi7s+o3Zryz9oBSgUHOyao5X0I9y6X7CoqW1/kfCpi0h4Np9XWOiu/93QYFFVsyxemeFf/WEDETKTlbX3K6OoaqD4amIwJLQufDDbt76eoK7EUG1RHWph/y02NDhRWYdMJShvGKjoK3wqkJ85hz6IrNhxXi7pi+q1o9a/4Cq4GiYCO8Xdo0Jrcbb1rC9IWR7kJST1lP+eyFkJo5OfKeLRVaP/9CYuj7WEdPoAJGHetsUOk21sZTEk91jown9JHc/r8kJGBca3GTlObUI4JylwMnSUokmuKruSgAXMgyn1ck1IKjqYOUuod0kv0vuuVK5hhSFWfpT6uxMkru4xHhXK8tPjAfKZPgyCR6o184CqiT4WAK8a9unZdtRye+uTx1gvSxC+tOQAN/4xU4ni5BDpFbl/RiBdUrbQbt7owl9geWtrzlYdqkzT2BZtMA4VJlNar9Q0DY0WGzc0gsRZlXosFxv93PH6syED02oseq7E3YmNFmvqbWXmB86WuI5eC33pzy2L8HnNXYFeFZ/J2R8+NdByEQM5hwMQczFakp413p/9FyYixUTWGV+lFQCyyczCfHO7XLHFFmVa6UzQX6poDcK+aUVVrYdf73OvNytzIUKnUOmas5Y2e2s3r/oCMgg9LmXh+US7X2aquADSeVlJfX7SMhKOSpcyFAhWT0DhwoHLUA6YKgwnmfkiYSwZENsRKGIGTEolaiprfcT3MXkTpn9rMhStiSDEuiNyanSGwEKmAR2ERQXNkTojTKU3IQM1daF/8z+stMQbZeFhuwKsyyq8KFs86peVm5Dh3nhrfPqP3nrXShxL2To9g3XAeW6TnWzJg4Zml0HrJcFMGRIVk8kZMhwISGroO9chdE2OvzZN1Q0rwksc+KwI61L52a5cJ7OvNDgRX3f4lJQbAT5lXOu7PJKBe/NfsVlcLxtz4QPq3Ck349D/d27LissGyeR7iKme7hZU0waTQgZD4osQkid8AHew1nTal+A6TKEWLlaZYiwzKmyocIL8+NCg6WQujR5VLsrQZEBuyuxocRAgF2qsg0zGlFq5ymFljfCUYfCSqn6iMNDgia8PxRAhBALi5ESkhJDJRF3cEvi8w02hBj9EYW1elNB8rsfJgSgbamFWogwU5W4ci6WTUYvLgQ6M++dE5VfCFRuEuQhZikadgQhsL0ryLYa+ZVAbQGdAZLbMKIoQNkQYRkWUwAKiBboLIMgh1aq8vvzop7w7nLOtHf9eW4KjboCpe4euFBZpGBpWZi0jaEKkzYVTCWETErvb3TmY5GzYo75CqemrcJ7l+t3uVbKE17eVDvaG73nfgCYBHVYgbWpXKxiUzlUbt/dNVBcALs7gN6YZXFh1vv75TbE6NrSGykT5Yss3g/tO29eEv/ePIut96Dl3iX++Q/y+8kRhuSM6KOD6GQRsnB6FRLt1XBLsrtrNyxx4BcBBSoHyyW3WxfLOVi7O1UpBiOwjHOlNwAKwe6uQHKguDRT4GjfaCvMj7o17tXmCQBl2sq22gopaw5BsHlaIN9kUFubIQ8Au8JMv6MyoMiBLDO1tEpHSdevz4nHLAOKAgIE0+3kR9znYATj3uYODhghJEkosggJWVjlau1XIwf2Q4WRCZwPlm0QMeGmcL0TIK70gd23HFHonyuzIwSduLIOV2GFV+HlYxUXgvzCJbzbU12Z0GF+ZUSUjplFBVBsALUD8jsmRKgvAbUT5BfeLdlpFBsxywuFbLuzZR3M3IjYFaa/ua5dh86UGWnoHK7C3oNYSE8X+4IpmM8wOpdhKKCCkOHiBNYB0UjIOdHL42WokJAZOXZkYdv2WMmB8ji1f07fySr0/gPVVXhXMOuzupjTGwUUusrB2njtC5BfmtSn/Eqwu2ME1u6uEVm7+xrFJVBcGSG2u6/N+rt2vzvmOG3b8QdVFxt7vsyIHb0JQnmZCxUG1wFUwtHPvQqFrKuSH+Lua1uZhAOfH0cYEpIeXfUQnSxCFszJocLQyfJrYzW148oX+LlYXthQu/Vl7Ss73+CFMonrG1UmuZuRhOa42wd2+ZxJaM+vgfxSo7gwTtXuQQF1I8jvFMieKhRXGlopE2K8BrJb43zpe4LLh2bd5SOTp5XdmtpZSgNaK2hdQCsTOiwuFGRXABcZZFeYxHrvmnWmTNJ8ZguPFvAS3COujYjNrdd7Tlbr58SQISGrg9mKhPjMGSo8Inl478EbCxUG26IFSGPsTZys6vdH6/1SB0qZAqImEcqEEbUpoSBaW0cJZU2rqvSCaSK/BiDA7q7pd35XQ19obJ/PoS8L5C/sgGuz1JeFWX+hkdv9d3fNXIa5TYgv2710oUh7ftsfvRHTv8z21/a/Fg4FylBniXOzave32/2sjdL0729t3/rneJTAmjMZfWEhd0LGovNfIUOF5CwI82kmPXfLkPshQoU2z6eqPB4+2IO27HbxXKpQaIgLGbofK8bEhgpFm+rs0KZcg2hTywraJLhDGyGUXwnyazNacHdXQytge7+A3mgU93LIvR2y+1u88NJjXNwxy+z+FnJvh+JeDr3RZn9ljt9dWyfsSozQqp3P9KO4sALQ9tP1txRP3nWV1+nwnKzyM+h4P6OV38cIGU5dwkFaRCMhK6SLLmK4kJAmYonfc1Ho/jk4DcnwnZysWD6Wa9NPeFf1tlyuU1kfy4YSdWZypYoLBa1scVEFbO+b8GJ+ZQXWPeNKFfd3QCF4x7sf4na3wSvPPUKhBV/14A389tvP48vuvYVff/QilGi8+vABLjc7fOnzzyF/QUM93tgkeSOqoAXFhYY8NI5UdquRXyqobWGn5zH9FBvilF1RLl2yu1YC0U5Eikkc23OfBCha6lPFkt+BSmwdkzCe0vyDrMtFyB6JPEEISZBUBBZwXJJzwwN7L/TUcp3R0g1mw97rKkRoxIIpHmoFlq2FJVqbeQi1DdkVJpEdQBnyKx7sAAGuXniGq4sdvvql1/A1z7+Gb3z5VwEA/84rnwQAfOPLv4qvef41fPVLr+HqYoerF54BYo/327tyoxDteS+qUKHrG8SIQDjnzYUOG661fO1Ps9N4EzuG/o5xf1JKfk/p74WQRKCTRYhPLSQ0sZPVdK6hQoVhTpbDv85Ym25iZFe6wd8vCBFqr2xDKaS0tqP6TFkGKcy0OFIAu2uzPr8LE8q7Z4TGneef4eXnHuMrHrwBAHjl8iG+8up1/J57n8G1bPG117+FL+yeg4LG490l/uUXvoDfyF7Ea0rj6VvXyO9qbB4bh2zzxJxn80QjvxBsbozQUzemX5kVXOoW5vp2hZloOrdu0y6P18pypRxi4sh3pdoqv4chw2OS35uS6qeu+q67Xwsh5wJFFiExUvqvvCFUuFcfq422aXb80KBf5yo83hca/tKfA9COrNNShQ5zN6pQmRIKxcZUa88vTSmG/NqUaLh91w5QRmD9S+94Cx96x2fx9fd+Hf/uvVfxufwW77u473UoB/AMn9n+c7znnZf4O2+/gl+4+ir8YvZe/BaAp8Ud3N4B1O4C2w2QPQO294Ds1tQclVwgGSCFRn6hkHmjDEUBKMRO3+PVtKrVAdOVk5VHCpA6gXEwhNh9Gp1GByyVkKH7m2HYkJCShJ4khJA+dBZY/RotH5YSiqlwP0c46s4VJ4WZykZczU1bUNTUyDLrClswNL/WQCG4eHCDV55/hA+947MAgP/gwZcAIBBYFW692+9D7/gsXnn+ES4e3ACFmHa987h5EbVC2QcpvCl3XFHVtutruBfV/VLx/U5klM+bEDIqdLIISZ0x8m66lm/Yc62C45yLZYuOal+gaMAVBM2vqv1MCQUgv5ayunt+pSF3cjx//xm++vnX8Puf+yS+5e4zAMBdddnaVbf9L777E/j7T67xcHeNxzdX+OKzC+TPFLJntnL8HQGeaajM9kOZ0YfZjVOBKK9Vi5hoGwSys9eUew6Nkqo6vLs/ZVX4Fmfp0PZjGKNNQsgg0MkiZEQ6T77bFp4c4wF6qM2yErxuXF/WjnLTzZTFSFHmZ0nhTWezgR1ZWE2pIztTcqG4LnD/hSd4z4OH+KbnfxkfvfP2UZf10Ttv45ue/2W858FD3H/hCYrrAru7GrKz3d+Ykg5auel4tE3A12UpB+0+MrHXZa+zVisrvH9N60Om/iw7hr2Tn8SckIVCkUUIacfP52oq6+Avgap0gy3jYEKFVoS5UGFmc7OujEh48e5TAMC/eec3T+quO961l1+ZqvGFPa87v9+vskq9f31t4dJaEVZ+jRJC4vDbgZAYS0zePeZh7zsdgZjYL/Vgw1KxJHi3i9ZllXezwvwYsQWv+CegN0BxoaEvC1y/9BQv3XmM3/fiZ/DezX1cSNb/WgBcSIb3bu7j9734Gbx05zGuX3oKfVmY82xgzmuLo2rlwpsow5qu+rvErt0tlarmbTx0v8rXR3w2SxRvS/y7IWREFvhXTMgEDDS6UHcNDx1R7X2PDqPUWguRBlPElDWyiiAXyYXOXGFSN7Gy9kJumRllWGyMiPFDc2qrkV/Z0OEdDWw0bp9t8MVn9wAAWx0ZrdcDd/wXn93D7bMNsNHI72ibgwWobT2EafpnR0XWqtRX1+XmaqyFSB125KH4hVrD++nfP3RMYu846nCIqu+df08PkdKoXEISgH8RhIzIpDlZHZyPxil1jm3Xc3R0ZD/JzXlUboqAaoXyBwD0xmzPLnJ88IVX8d7LN452sRwXkuG9l2/ggy+8iuwir53HP39xIVB5vZ8+tWKkrcnsPb5Gw6l12ujaLnOyCEkWiixCRmRSJ6sDnSeH9nE1okLKKWaM0JJwv1yb6WlyM2+g5ChDh6ZuFqAzDWwFWabxxu1dvCt7eNL1Od6VPcQbt3eR2fZ1Zpws//x+v/RGmdGCDv96/OsMCZ2+Q/Rxsnq2GWVqJ4sQUoMii5ARSW10YauT1TgJsaqHyPwJk91Sa+P85CaUKIUVWMHoQsAUIFW3AATInirIvR2UKvB1z/0m7qmbQa7znrrB1z33m1CqMJNLPzUFUdUtyrIRe6MLbX+1Lc+g/YmiC12/bsALlTZ8drH72cfJ6gqdLEKShSKLkNQZw8mK0VhR3LohYfmGsG2tgcyVQ3DFSJ0wA5Qto6BuTW4UAOR3izLpXMmwSdNle2LPA5uTdWu371Ce2/VTZ0ZoweZmhUTLODQ5WQemyBkMCiRCkoUii5DUGcPJitHmZAH7TlbYtkjlZNnCncqF4DRqTlZmDavsSTUSsdDDfh2V7Wl7Hpjz+k6WO7frp+RF5WRF7seekwX0c7IskzlZhJBZocgiZEQWlZPV5mT5yd9hIVI/J8s5WUogOyNYylGFnpNVXALQQH6ngH57g6JQ+MTDL8fbxdUg1/l2cYVPPPxyFIWCfnuD/E5hhN5l3clyoUy/v6WTFcnJ2nOy2nKyTh1d2BXmZBGSLBRZhIzIonKymnC5SSGxnCx/v8wKlwzIttoUAdU2rKjNnIGSC3ChkeeCFy+f4Av5cyddn+ML+XN48fIJctu+5GLmUPTO7/dLdkYgluh4TtYebTlZMZiTRchZQZFFyIgsqk7WMe16zo5E9nNzGhaZqU8lBcofAGZeQAD5NsMn33wFn719cZA6WZ+9fRGffPMV5Nusdh7//GqrvSrwsdCgm9NQt9+vsUYXLrFOFiGkBkUWITEGqly9qDpZ/oO/0Gbf0KnxR9lp7YXMzPFic5kk10BhRu9B10NzxYUguzFCJ3sqwE5web3DO6/NfIVD1MkCgHdev43L6x2wE2RPjZOV3dj6WH4IU9t+Fqj337su2Outjap0WMdLu/sRu5/+/cOK62Sx4jshNSiyCImxxMrVfRwVh/9QDB78e25LORl0Ue0bHGMESjViEBpliE4K2PpUZhJm2QFqK5BbhWev38HrT+/j5954Hz67e3y0m7XVOT67e4yfe+N9eP3pfTx7/Q7kVpnz7GDOm2vTj6IKHfrTAEle7Ce9+9dbFNG6WdH7Vb4+4rM55vOcmyX+3RAyIvyLIIS04z/sY65JRHA5p0us+yM5yhGJYvWTygG1BbIbI07eeHIHAPCzT7/8pO6641172Y1AbeFVd6/66Pol2nOqWq6rxN9viWKIEDIJFFmEjEhqows7t+mVa2haXwuduTn+ADMJswsdBqMLpQCyWxueA6A3wOaJQD1TePzmXXzu0XP4ybc+gB97eu+oy/qxp/fwk299AJ979Bwev3kX6pnC5olA2/IRaqeR3dqcrHB0oQsVuo/MzV3YVIjUp2uu29SfJXOyCJmVzdwdIIQcoDiQeD1Gm7aMQTm6zh9l5yeJu8mRnUBRuso3EpRhuOymQGGnsJHc5kI91VB3xIzyuxFsn2Z46/E1Pn31Mj6++SA+/hD4i+/+BJ4Ut7irLhu76rb/mc9/HQDg02+9jLceX0M/zbC5cecEsqcayr6W3LhWaqtroUITOrTOVhlO1PtuVhG4XF1HbY4hZiiQCEkWOlmELJRBywBUjZbuh24Ll/nrfLfLhd3KYqRFORm0Cc2Zwp9lyHBrltkzAZTG9tEVXn3rAX7xS+8FAHzfo3cAAD6zfRztrlvv9vvFL70Xr771ANtHV4DSpl3vPJLb8xdV2FAr00+zvdgXVeH1NdyL6n4V8f1OZJTPmxAyKnSyCInhwiwpJPI2OE69CloWRfNoNV0AhTLnaWrTHS+Bu+W2IYNobVYV5rbJtkB+oZBtCxQXCpIbJ0nvbI7UtprmZncXuPzCBrv7Gk/VNX4LwDbP8OrNA/zwG/8q3nf3dXzl1et41+YhrmWLZ/oCX9g9h39x87X4zJOXAAC/8ehFvPbwPp6+dQ31dobNY4G6ATZPTDcv3ra3c+eEnul/trWCcOumD0LlZpXXh33XyiXBx3D3p6hEa+N97YiIxIWWkjTcLI4sJGQPiixCfJyAAKYXWLqIn7MhtKe1PjwXYZkrpLzVwXH+Of1jyu1SF1gIhJZSQFFAZ5kRKBew7xVkp1FcKKjclG5QW43dtSkCursSbJ5p7O6akg67+0D2RFA8AJ6+dY3PF4KbPMN77j3Eq7fPoYDgp978Gnzt/d/GP3v8ZfiKO2/gC7cP8Hh7hc+9/Ry+9Ogubh6ZivHZEwHETt+jgM0TDa2AzY1Zqq3Jt1JbE+ZU28KOPCyMv791eWV5vaq7L7DcfSkizlXNCVM1AVITSu748Jjox9mwvUlgTS16/Ouk60YIAIYLCWkmpf/Mj3EqGgTYfqmB5ussa2WFQiLyWgojUFyNKXH5TdY1kp0TNKaEg0k4B1Q5j6EN7T3aABq4efMaN9sNPv36y/jUWy/jp197PwDgx1/9IADgp197Pz711sv49Osv42a7wc2b16bm1aNNvb0bGxLc2fNaYSW7qm/GuUKZXwYVFCNtum57X1pDecH9bRTGx0y1k4KD5Ujp74WQRKCTRUgTKYQKHcckvvuulBMMSlVOVlvyux9CDNt0x9mSBzqrEuJlV0BfKrPcKDuPoRUzmSC7LVBcZNg808gvBRePgfzK9lMAwE5/8/gC+bXGm7cvQF8V+LVHV1AXBd58cgdvP7nCa3fv49Hb1yi2Cvo2w5MbheyZlDlY2VPB5plxsrJbjezGhAmzWxsivHU5WF7pBq3N9DpAuXQipiy+6tY1lXVoDQ0GuVthqPAY9yeVUCGw59gRQgDpmEypAeAj6tvG7Q0hc7MXKjtdaA1S9R3olpsVy7sKQobRqXXcuf28rKyqvC4i5r3Y7fZYfeFqIyggU0bEZWJChwCK641dmvf5pYLeCPIrc77b+8quB/IrQX5t8rTyaw0pBLcvmId2cVkAGw1kuhoFmJtq8erWtHH5poK2ye7ZDZA9A7IbU7IBAC4fm7aymwKy06XQUs9yuzQzWEue24r1BZAXpRiSrZ3h2oksrYE8r7tYucum194IxP06Y3tCKyraGkKMtX0OfH9PWcIhPBdDhmTl/MPi+93Lxi/5hP5VJyQB5nwwHOEC7IWeQnekNlVOMGqwT7kBl5cUjLqrFe8srDAp/NChqf5ehgp3ZmShKh0ltzRNZM8AaFM7CzAhP9kKLt7KILcK2Zsb4JlZyq0y67dShgY3TwTQth2/XXsedWvPv6tChq7KexkiLOy98q6tViPL3MSyfEX0frXcz6iTdeBz7DXIoezjjK4SBRYhABguJGTR9Ep+j1GG/mzSvSssGpsixq2zNaQAXauRZZBylKFsTdK4bAsUIlC3uQkhSgERBbVFmSNVXAguHwH5BXCVA8Wl2NpWYkYg7gApFNSt4OKRAjRQXGps3hYUG1RzId4YIaW2QLY14kptrcDautywolpqO6pQa9NfoBpVGJSjQB5xnZpKOoQu1kHH6cikd0JI0vRysjxrjBAyNQ0P6oMP4LbtZQX3yFeBLvbP6btZbv4+X5A5IeImVs7rrpnszLx/roSCq/xuzmecJ9EmzLd5akTT5ompabV5LFC3JpFd3Zr3kpvt2Q3s/iaBPbsF4DWtbJK7yo2oDHOu4E8I7V8H4I2i9EVm4BI15WP5Ydgmhh5VSAgZna56iOFCQkIW5hoMEjKM5dPEhJdfysGrel4W8fTrZ9lReigKSF6UITdXj8oJLLUzblO2NbWplJ32JnuqayIqu0FZX8vMeVgXYdlTU2RU3Ro3Ktuadv3zANa5sqFOcTlXCpUg9K6jvC53r/zSDXsfhIqUcqiLuclChXOysL8fQsaE4UJCFs7JIcPG46oQolY2szPPqwR4vzJ6KbyMUyV5boTINofOMohowLlHUMiKHEWmsHlaoMjMiELzY/K28guxQgkoMqCw8x4WF27OQSO0TM0rM/mzFCZsmW2t2LICS3bGwTJLmzOWFyZcWGj7PrdT6ThR5Ikrd53+e5vgrpsS3HvdZ4YKCVkrvZ0shgzJWTFA8nDyk++2iYQu1+9ylwrP8fEcLVdvSgpd/gDVVDalAPKcJ5eo7vbdPDNiavMUkJ1Zqq1Z7++X3erSGVPbSmD55wv7UauH5UpbuDZ9cdV6D1ruXeKf/ygjCwlZMX10EJ0sQlKiqep7X9qm0XGnijlg2kyRY3fw6mxZN8uNMnRtezWzymO0BnY26R2AVgqC3CxFwyQ+meOVckJEoLUySfNKIBuUtbQ2TzWKjamjpf2qExqABi5sUdPs1hyrdlZE5SYHzIksseFCyU3iO3JtBFa5bHGvanWyvJCrL1BqCfH7oqOTI9Vjmp1WKHoISQLmZBFC6oTOxqEHdljqIKj9JL6jVdiRfc5V8vOldOVAuRpWmVfqQXKNzY2Z/3BzY95ne6UgirIN6Cofa09gFUXlsPnCCrDbvHvQtxZV4s4VIWQ6KLIIWRLHjjCstdHB5Yidp9D1aXbMicOOtC7dpMxOaEleQMoJmuv7Khv6y24CIXUTCLAbGyK8DY637cm2qM6V1/txqL971+VPo9Nwjw7Sw63iyEJClk2viu8+rP5OVs8IE0UPUv29S+V3YD9c2FAZfq8CfKzsgMrKdeKqvvvtbLLqHCLQmapG4SkpK8G7fbRSQCbQyu4LQG/MsriwbQpM9Xixx4jJ3yoyKZfl/IjaCif7TeVGMJbT5LjRjS406Is+V3jUri9HFLrvxp2t4u4Jy5rIKvLqPjXUx6p9z/oiK/z+DQTYUSJryirv4TmZpE9WTkM+Fiu+EzI3nQXWFMSmbGlKgBfZT+IOKr+H68ukcU/AlO6RI9dVGQXvGJUXppSDLRqa3ZqRgOU0OHkVGnTrTc0te5wVSWViuyewapcVE1jh9cUqvYf3qmnC7UMCa2aS+n0kZKVQZBFyiIGSiAdzDrqGDI95oMfaDtrVYUitDCFGhJbfjyJIOLeIhhVUNl9r64UQd0FNra0VVHu1tqxQs8e6cKSyVd3Lc7lEdyewvP5FBVZYugGR+xxzb475rIdwsXow2O8jk+wJaeTocCHAkCFZOb5DMVDIsJd7MEfI0G87Nmk0UE4cXR7jwoLB5NFunXbH+6FDd6wNEWoRIJMyLAgXRvTrcCEyuhCo7WOKi+oyjIhcV4ntTogFDlvpYNXqZHkCywsvlt+XuRcijNXK8gRMiqFCs+sIIovhQrJiWko3MFxISG9SfmB0fUCGblbDNfWamscPg4WhNCc4vPU1R8uFDoOin2WeVFGU+VWyza34QSl6lC0gqnZFdS5ta19tc3Ocjrdr+lLUBZbfP6AKfYb1sUT2w6ot9/PgfT3WdUw54T3lvxdCZuIkkcXCpGTV1ATECkMifsju0CTSQFRo7eVmAVXdrEiozU8eFxc+rCWho56AXni5Vtu6KFLbKncLhZewXriq7tW5JUh43xsd6fdXeQ6do01gtQif2n1NKB9rMPy/i6VN/0NID47VOyxGSkgToVOzVNoKk9pttcKkZYHRhsKowXQ70LryypUyRT1Fav/CSQGzL1BNvmxDi2JjgHqjIAVMCNEmlEueQ4uZdgcFzL7uc8l1uV6VYsmee1cAqhpdWLlekbwxX2CFYUIA2g8hunPEaAsVtgmstYgvOlmE7EGRRUgbcwqtI6q/H5zH0OxUv6Y+Qssd660X8YRPUZicLa39wu5mvwImP6tsQ8PUadBWUNkwpBVI5Xqtq2SsXZmIZRY59t0yl4fl0qZiifpoEVjunth9xQ8VNpUq6COwOoiRo+YrnNNtpcAiJMpJie8OJsCTVTNwAvxgye9ANAG+UWQdWzfL74cfSvPW1c4ZJsK744LaWmVCfJgw795rbabi0bpKlPcvwRlidrsURSXeYiME/RINReBwuddBortZresuVhgm9MXNiXWxaucM6Vt5vnXXAUQRE97JmdAhVDhu4jtzswjpTq8H3BHuxFEuSIy2iuZa19yb6DmbXCLfSQJsDlVhwoiFtksvUd0uRWuIrZ0lNk9Lgv2i7TSVaAjDg5FrOCiwDt2vIxjbxUp+wnJCEuJUfcPRhYScC6fUzerwEK/Vz2o6VyBoXCHSvTBdXlQ/QT2r6NLf398WO0+DoCqP0xE3qvGi98OEnVlLLhYhpJFBwoUOhg3Jahl4ip2xQ4ZAQ9iwLWQYbO803U54X8LQYZYF+0XCgkFbbpqdaB/9fK7Y6D93SCje3OswjOjv49W+2nOwwrb6Tp8T6WNMZB1dgHTO+lgMFZKV0sPFYp0sQgYhtervQPcK8MBRD/62QpuxfaNCIwzR+cU/wzBiOHdg6D61rdd6X2D55/HP7e8TTjN04BoP3o8pBVYPWOWdkGnh6EJCUueIUYbAgdpX1U7H9aUcfeicpXofy5GKfujQJcOXIxG9NsJuOZ3k5oouNLSSxqXbx548vozmZNUF2f6UOWE4s3v4tH5B8eus7yLNQqsNCh5CkmXQcCHAkCFZKSNMsWOaGneaHSAitmI1s7pMt+OfJzbaMOynEz8S7OPa31sfCUM2vffZE0axsF6DcxWKq5iIatoWuEKdRxSG28Nj9/Zt2TZXwjtHFpKV0zPhfbpwIUcaklUywoNk0LysPg/QQwnXYegsRig0WkoO7E0oHTt/zG1qEk/h0t/e5FrV+nVAYLX2q0PJhR739yCHPtcegr/X71tXKLDIChlSxzAni5Cu6PYHbe/mRi7lUB7aJTerZXtvodXWV19ohTlVudeenztVy6lqWPr5VkXQng7a8fvRdl0TCayTSm4w4Z2QpGFOFiEzogvd3WE4lJvlcpyO6oh9YDaE5cocq/AcbTlaAFC4wqPVIeXRoeBQyjtel1PnlP0Kt/lL/zr8dQ0J6HVh1CNEGBzTKpBOFSEpjigkhPRiFCeLIUNCFkBb2K71uIbk8Fiory1/KXSyYiMR99rrsa3Yd7L2RFGsTENHgdXKMe4hIWR2htYvgye++zAJnqyO1BPggdOT4M2OjfvU2ohNveOvb0mIb+xTW5K735c20aIjos6ni3sVbo/NSxjrS2MphoGS3cN+HYAJ74R04wSBNU+dLDpaZHXMnQAPHH7AnpoEH4qckyrFN7haXu6UDh2m0OFy6wBTMLQoqsKhB46rtR3mbLW5V72uNbg/MZE4ZLJ7z/w8JrwTcpix9AoT3wk5lhQLkx46V5cw1oEHaDSfCegWVtNF3B2yRUxd27URiW2hwmAf//haUdGwn7HPrikMGmw76P63uVxd2xgQFiAlZD5GDRc6GDYkq2OEsGFvx2HIsCFweMqdYJ/GGlqx/jXVwIpdg4ok2HfFHdck8MrXxwnFgyMJTxVYA7tYowgsulhkZQzgYs07rQ7DhoQcZhGOVlDaodHVAvadoJqY8UJ2TcLGX4bHx9bHXKvwHGF/Owqs/WsdQWANDEcUEnKYsfUJw4WEnCGNU+70DB2aXToKrbBNHQgc/6epvS7rmtqKnbepn33Cg7E2G3KwOk11RAhZDZOECx0MG5JVUU4LM+z/KlOFDVsf+D1Dh9H22sKHTfsNKULaxF+5T4/wIHC0g9XY3qH+lQfOFCYMz89QIVkRA7pY84YLHQwbklWSekJwwwO38YEfK+vQse7TnqvVFo6L9TEcIRh1nVrW7dXeahEFLeHMzlXyY/s1lMU4WmDNTeq/34QcwVR6hBXfCTljykruPkXRLLRa9m0VEc6t8h/Yztnyc6lq1eSPFFr+ef1leP5wv73meuZgtew7ZS4WISQdJg0XOhg2JKtihLDh4CFD4PTRhtXOnfc9GEKs7dwy0vAQXUYUhhwKDZb7jSiwOlWPnzFUyDAhWSEjuFhphAsdDBuS1eCLiAHDKr0flF3O3dJmr0mkw0Khbfv26EMfl6nTfj0EVvN+kYEA5yawgGFz5QiZkan1xyxOloOOFlkFKUy10+f8B9ptTIiPuVUdXS0RiYcmO/RnECLiw/XnpPyrBnF58Ht1BIFlDhlJZNHJIitgRIGVlpPloKNFyAwM+iDulxAfFzS6/jMEB9ps78+RE2f36RshZFLm0husk0XIqYT1noZq9piHcdfz9w0dAu3hw9i+LS5PJwe9qbjosfsdOndTn3vej5NDhABdLEJWwqzhQh+GDsniOeewYXVQy/mak+MbQ4kn0hoSBNpzyVrFEsOEhCyBiRysxi8vlnAgZChiJQ4WTKPw8QVGKJza7oE7LjgmDN0NIbZ8sdMrqb86qPcxqyvTsLbrIWQGkgkXMj+LkDijhQ3Hyg066OYcFivhz5D7HyWwhmAJLhYhKyIFXZFMuNDBsCFZPCNNt2OaHKF+FnB86BBoDx9WDRw4v2ougnosrr0hRNUxU+WUx46ThzWKuGJdLLISJhZYaY4ujJGC8iRkEFKYjmTIRPxj8pqqg9u3uzb8trrW3ortH2vvmH61tDFoeHBlvyuEzElKOiI5J8tBR4ssmpGS4E1zIyXCO05xtYDDblQKeWsDhDTbj+/xlZlCiJDJ7mQlzCSwluNkOVJSooSkxKilHbo2d2jk3aFRe23zEo5BOIF0Ewf6PnhyewoCi5CVkKJuSFZkAWneMEI6kaIbMEANLeBAIc9aOx1DdU74DC5gdD9B19G9misHaxJS/L0lpAOp6oVkw4U+DB2SRRKGxVIIG/bpR4/2O5ddGDKx/VQ65nz1cq9GFlijhgkBiiyySBIQWMsLF/okcAMJ6U/4wErFuRjI0TqKvsnsYzFGP5bmYFFgkRWQuj5YhJPloKNFFseICfCmyROSyAdMhi+b7JPUPoer1UNcdXawRk5yrw4dsWQDQJFFFkdCAqvxi29RIgug0CILZOSwoWly5NAhME74cO8cA9TL6lofq4FRwoNAOiFC02jwPpmveEI6kZDAAtYksgAKLbJA1uJoAdOIrfJcal80xV6fIKoco4krgA4WIQOSmMAClp6TFZLgDSZkVk56EI+YI3RyyYNDRUq7Fh09wKjzDqYmsAhZMEt7/i/SyfKhq0UWRaqOligjBkZ0tcpTJVCM9ChRdYx75e7rEdDBIqQicXG1LieLEDIw5Zx1faey6f+g7jSJ80gcfe5jw4OpjCQkhMzC4p0sBx0tshhGdrNMsxPmaDlOOWesG0c4XoOLt2PdpBRDhHSxyAJJ3MFyrCvxvQkKLbIY1iq0gMHF1izMlONGgUVIxUIEFtAisjZT9oIQEuGYXKixcQ/lY/rlhMISxVaiAwhOItV+EXIGrMrJ8qGrRZJnAjerav5EwXNq/1IWXKe6RyeKmNFHENLFIgtiQe6VDxPfCUkO/4G3dreh0ONM03MKKfZpaCiwCJmV1TpZAN0sshAmcrROdrOqhoZpZw53ayhRNZAoHtXFosAiC2OhLhZwLonvTVBskaSJjaI7J7EV4vpX6H5CzN9/igTyk5qZuPgpRRZJmAWLK8d5iyyAQoskzgTzG1ZND+ggTZ2w7/d96lDfgCHdSUUWBRZJmBUILIA5Wav5IMlamfBBOOgDfupcMpdHRYHV8WQUWCRdzuG5fDZOlg9dLZIsE444NKcYIS8qtXIUxzCCeJxEXDEPiyyAFYorOlk+K/yAyVqYeMShLvT4c+QtiZHEFQUWIYZze/6epZMVQmeLJIlztSZyhkZxtUzD47Q7JCOJwslCg+Vciav+qiYL5QyEFZ0sQhbLRK7QqFO6pOpsjdi3yQUWISQ56GR50NEiyTFxjpY5zQT1q+Z0tyYKw04CQ4QkYc7AwXKwhEMfKLZIUswgtMypZigWKqo+l2PXeR1jx8zg8Ew7cpACi6TJGYkrB0VWXyi0SFKck9CKEYovYFYxFYMCi5CzFFgARdZxUGiRpJhJaJnTJSK2EmRScWVO6L0+y69mkihnKrAAiqzToNgiSTGj2DKnpOCaXFiZk3qvz/ormSTGGYsrB0cXngJ/gUhSTFxLa+/0cwiMhKDAIqSCz8d26GT1hK4WSYaZHS1z2vNxtWYTlxRYJEEormowXDg0FFskCSacWPoQaxRcs7p2nOiZJAjFVRSGC4eGv2gkCcIH74wj7dYWRqTAIqQOn3v9oZM1EHS2yOwk5GpVXViOu5WESKS4IolBYdUJhgungEKLzE6CQsuRouBKQlg5KLBIYlBgdYYiayootMishCILSEpo+cwhupISVT6xMC9FFpkRCqxeUGTNAQUXmZUERh8eyykCLFkh1QRHD5JEoLA6GoqsOaHYIrOwIFfrLKF7RRKB4upkOLpwTvgLTGYh9sBOZJ6/s4cCiyQCn0/jQidrYuhqkdlIOCn+bGByO0kEiqtBYbgwNSi2yGxQbE0PxRVJBIqrUWC4MDX4i06SgSHEceH9JYnA58700MlKBDpbZHKYGD8uzLsiCUBhNQkMFy4Fii0yKTGhBVBsnUKTc0WBRSaE4mpSKLKWCAUXmRzmax0P867IzFBYzQZF1pKh2CKTwjBiPxgWJDNDcTU7THxfMvwDIpPSJBCYwF2HYUGSAHw+pA2drIVCd4tMgogRDczdahdV7j4RMjIUVUlCJ2tt8A+NTIITDhQQcXh/yITwe3950MlaAXS1yGQ0OVrAOl2tthAphRWZCIqr5GHi+7lAwUUmZW2ii6KKJAKF1aKgyDpHKLjIZLSJLSBtwXUooZ/iikwEhdViocg6Vyi0yKQsSWxRXJGEoMBaNBRZxEDRRSbnkPDa238AIda33AQFFZkYiqpVQZFF6lBskcnpK7amgOKKTAzF1SqhyCLtUHSR2RlDhFFEkZmhqDoLKLJINyi2SJI4AeYX/vTXEZIYFFdnBUUW6QfFFiGE9Ifi6iyhyCKnQdFFCCH7UFQRUGSRIaHgIoScMxRWJIAiiwwPxRYh5JyguCINUGSR8aHoIoSsCYoq0hGKLDItFFyEkCVCYUWOgCKLzA+FFyEkJSioyEBQZJF0oNgihMwJxRUZGIoskiYUXISQKaCwIiNCkUWWA4UXIeQUKKjIxFBkkWVCwUUI6QKFFZkRiiyyfCi4CCE+FFYkESiyyHqh+CJk3VBMkcShyCLnAQUXIeuAwoosCIoscp5QdBGyDCiqyIKhyCLEh+KLkHmgmCIrhCKLkDYouggZB4oqcgZQZBFyChRhhMShiCKEIouQQaHoIucKRRUhe1BkETImFF1krVBUEXIQiixCUoBijKQCxRMhg0GRRUiqUHiRsaGgImRUKLIIWRoUX6QvFFOEzAJFFiFrh6JsfVA0EbIIKLIIOXcowtKDIoqQVUCRRQjpBsXY6VA8EXJWUGQRQgghhIxAo8janNoAIYQQQgjZR83dAUIIIYSQNUKRRQghhBAyAhRZhBBCCCEjQJFFCCGEEDICFFmEEEIIISNAkUUIIYQQMgIUWYQQQgghI0CRRQghhBAyAhRZhBBCCCEjQJFFCCGEEDICFFmEEEIIISNAkUUIIYQQMgIUWYQQQgghI0CRRQghhBAyAhRZhBBCCCEjQJFFCCGEEDICFFmEEEIIISNAkUUIIYQQMgKbuTtAzoOX5N36FrfmjQjEbdh/AYj/OrY9eB/uF2yKrmttA9DRdvb3i7Wvo/s0HQtAxBzTuk/Qfsd9a/v36RNarqNDOwf7eOjcfc7buE33/ByCY7310WbELfT+pqCt2ttaezqyf9Vm+GtaW19u1972YL/y2PC4/fPW2676LHvt6712atu8++JfVlP7/jma2q9db6wNAdy3Srivee1vq75/fuGXbn5Ua/0HQchIUGSRSbjFLX6v+ghECSCqXEK5b2SBKGW+9culQMSarUrMOn+7VMdW21V932C7FjH+bbA9XF++97eLfe32gxUDdps7ptzX226W5jrK15Httfew61RkW+34/aXZJvV1e9vRvh376/3tXdqOrm9r1z5so8eH22t90NU+8PYNt7uPrra/eV2KCfdaqn1FPPEg1f5i3/vbVfle2189XbatUF/nL6vt5n3tx4qN2Hrzuqi2e+uy8n21PbPvlWi73d/m1hdQsEt7rDuuWhblMQDqx8Cc2+zjji/K82XueHtdmWu33KbrbZfvXf/tOvtxZwJkECgAmQgUxL4Xu12goMx6Ma8AIHvPr7wEQkaE4UJCCCGEkBGgyCKEEEIIGQGKLEIIIYSQEaDIIoQQQggZAYosQgghhJARoMgihBBCCBkBiixCCCGEkBGgyCKEEEIIGQGKLEIIIYSQEaDIIoQQQggZAdF6f94tQoZGRP4pgGdz94MQQjyutdb/ytydIOuFcxeSqXimtf7w3J0ghBCHiPz83H0g64bhQkIIIYSQEaDIIoQQQggZAYosMhX/+9wdIISQAH4vkVFh4jshhBBCyAjQySKEEEIIGQGKLEIIIYSQEaDIIqMihu8SkT8gIv/53P0hhJwvIvI9InJfRP6UiFzxu4mMDUUWGZtvAfCG1vonANyIyL81d4cIIWfLtwP4eQCfBvAR8LuJjAxFFhmbbwDw2/b15+17QgiZgz+mtf6A1vqHwe8mMgGs+E7G5p0AntrXTwC8OGNfCCHnze8VkUsA7we/m8gE0MkiY6MBZPZ1Bv7OEUJmQmv932it/wGA5wG8D/xuIiPDXyoyNl8EcM++vg/gjRn7Qgg5U0TkPxKR77Bvvwjgg+B3ExkZhgvJ2PwkgC+zr98D4Gdn7Ash5Hz5FIBfsa9/B4DvAb+byMjQySKjorX+UQCXIvKHAFzZkTyEEDI1Pw/g20Xk2wH8htb6vwe/m8jIcFodQgghhJARoJNFCCGEEDICFFmEEEIIISNAkUUIIYQQMgIUWYQQQgghI0CRRQghhBAyAhRZhBBCCCEjQJFFCCGEEDIC/z9+Rkt7V1t5dQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 612x388.8 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"hp.mollview(template, max=50)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment