Skip to content

Instantly share code, notes, and snippets.

@smsharma

smsharma/TheorydNdJ.ipynb

Last active March 5, 2019 03:22
Show Gist options
  • Save smsharma/7db7271a4d36916ae6748489bdc31934 to your computer and use it in GitHub Desktop.
Save smsharma/7db7271a4d36916ae6748489bdc31934 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
TheorydNdJ.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The autoreload extension is already loaded. To reload it, use:\n",
" %reload_ext autoreload\n"
]
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from skmonaco import mcquad, mcimport, mcmiser\n",
"from tqdm import *\n",
"\n",
"from units import *\n",
"\n",
"%matplotlib inline\n",
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 1. Cumulative SCD calculation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The cumulative SCD can be calculated by integrating over the mass and spatial pdf of subhalos, with a constraint on the integration phase-space:\n",
"\n",
"$$N(>J) = \\int_{J(M,\\mathbf R) > J} dM\\,d^3R \\frac{dN}{d\\mathbf R dM}.$$\n",
"\n",
"In practice, this is done by integrating over the subhalo mass $M$, the distance from the GC $r$ and the angle between the Galactic plane and line-of-sight from the GC $\\theta$ convolved with the subhalo mass function $\\rho_M(M)$ and spatial volume distribution $\\rho_V(r)$:\n",
"$$N(>J) = \\int_{J(M,l,\\theta) > J} dM\\,dr\\,d\\theta\\,\\sin(\\theta)\\,\\rho_M(M)\\cdot 2\\pi r^2 \\rho_V(r),$$\n",
"where $l^2 = r^2 +R_\\odot^2 + 2r R_\\odot\\,cos(\\theta)$ and $R_\\odot$ is the Earth-Sun distance. \n",
"\n",
"Note that the integral can be performed either in Galactocentric coordinates (with $r$ as the integration variable, as we do here) or Galactic coordinates (with line-of-sight distance $l$ as the integration variable). Properly normalized, the answer is the same in either case."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Specify population parameters:"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"R_sun = 8*kpc\n",
"\n",
"# Minimum and maximum subhalo mass\n",
"M_min = 1e4*M_s\n",
"M_max = 1e10*M_s\n",
"\n",
"# Minimum and maximum distance from the GC\n",
"R_min = 0.1*kpc\n",
"R_max = 260*kpc\n",
"\n",
"# Calibrate number of subhalos\n",
"M_min_calib = 1e8*M_s\n",
"M_max_calib = 1e10*M_s\n",
"N_calib = 150."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"def rho_M(M, beta=-1.9):\n",
" \"\"\" Subhalo mass function\n",
" \"\"\"\n",
" A0 = 2e8/M_s # Arbitrary, to prevent large numbers and overflow\n",
" return A0*(M/M_s)**beta\n",
"\n",
"def rho_ein_EAQ(r):\n",
" \"\"\" Aquarius subhalo radial spatial density, \n",
" Einasto profile parameterizaton\n",
" \"\"\" \n",
" r_s = 199.*kpc\n",
" alpha_E = 0.678\n",
" return r**2*np.exp(-2/alpha_E*((r/r_s)**alpha_E-1)) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"J-factor calculation:"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"def J(M200, c200, d):\n",
" \"\"\" Annihilation J-factor for NFW halo\n",
" \"\"\"\n",
" return M200*c200**3*rho_c*200/(9*d**2)*fC(c200)\n",
"\n",
"def fC(c):\n",
" \"\"\" Helper function for NFW annihilation J-factor\n",
" \"\"\"\n",
" return (1 - 1/(1+c)**3)*(np.log(1+c)-c/(1+c))**-2\n",
" \n",
"def c200_SC(M200):\n",
" \"\"\" Concentration-mass relation according to Eq. 1 of Sanchez-Conde & Prada 2014 (1312.1729)\n",
" \"\"\"\n",
" x=np.log(M200/(M_s/h))\n",
" pars=[37.5153, -1.5093, 1.636e-2, 3.66e-4, -2.89237e-5, 5.32e-7][::-1]\n",
" return np.polyval(pars, x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, calculate normalization pre-factor to give the specified number of calibration suhalos:"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"def integrand_norm(x):\n",
" \"\"\" Integrand for calculating overall normalization of \n",
" joint mass-radial distribution pdf\n",
" \"\"\" \n",
" M200, r = np.exp(x[0])*M_s, np.exp(x[1])*kpc\n",
" return 4*np.pi*M200*r*rho_M(M200)*rho_ein_EAQ(r)\n",
"\n",
"norm, norm_err = mcquad(integrand_norm,npoints=1e7,xl=[np.log(M_min_calib/M_s),np.log(R_min/kpc)],xu=[np.log(M_max_calib/M_s),np.log(R_max/kpc)],nprocs=5)\n",
"pref = N_calib/norm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now specify the cumulative SCD integral. The integration is performed in log-space for mass $M_{200}$ and the los distance $l$ for better accuracy, and linear space for angle from GC direction $\\theta$."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"def rho_M_r(M200, r):\n",
" \"\"\" Joint M-R distribution pdf\n",
" \"\"\"\n",
" return rho_M(M200)*rho_ein_EAQ(r)\n",
"\n",
"def integrand_Ncumul(x, J_fid):\n",
" \"\"\" Integrand for calculating cumulative SCD by integrating \n",
" over the mass function, los distance and angle from GC direction\n",
" \"\"\"\n",
" M200, r, theta = np.exp(x[0])*M_s, np.exp(x[1])*kpc, x[2]\n",
" \n",
" # If J-factor is greater than specified value, add the corresponding \"dN\"\n",
" # (impose phase-space constrain)\n",
" l = np.sqrt(r**2 + R_sun**2 + 2*r*R_sun*np.cos(theta))\n",
" if J_fid < J(M200, c200_SC(M200), l)/(GeV**2*Centimeter**-5):\n",
" return 2*np.pi*np.sin(theta)*pref*r*M200*rho_M_r(M200, r)\n",
" else:\n",
" return 0."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Two ways of doing integral:"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"# J-factor values to calculated cumulative SCD for\n",
"J_ary = np.logspace(16,22,30)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### MC integration"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"# Number of points to use for MC integration (more points = better accuracy + slower)\n",
"npoints = 5e5"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"def Ncumul_MC(J_fid):\n",
" \"\"\" Perform cumulative SCD integral with scikit-monaco\n",
" \"\"\"\n",
" dNdJ_val, dNdJ_err = mcquad(integrand_Ncumul, npoints=npoints, xl=[np.log(M_min/M_s),np.log(R_min/kpc),0], xu=[np.log(M_max/M_s), np.log(R_max/kpc), np.pi], nprocs=5, args=[J_fid])\n",
" return dNdJ_val, dNdJ_err"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "18312c5aa997490abddea0460f1cacfd",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HBox(children=(IntProgress(value=0, max=30), HTML(value='')))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Ncumul_ary = np.array([Ncumul_MC(J) for J in tqdm_notebook(J_ary)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Manual integration"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [],
"source": [
"# Number of points to use for manual integration over each \n",
"# dimension (more points = better accuracy + slower)\n",
"nlogl = 30\n",
"ntheta = 30\n",
"nlogM = 30"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
"def Ncumul_manual(J_fid):\n",
" \"\"\" Perform cumulative SCD integral manually\n",
" \"\"\"\n",
" logl_integ_ary = np.linspace(np.log(R_min/kpc), np.log(R_max/kpc), nlogl)\n",
" theta_integ_ary = np.linspace(0, np.pi, ntheta)\n",
" logM_integ_ary = np.linspace(np.log(M_min/M_s), np.log(M_max/M_s), nlogM)\n",
"\n",
" measure = (logl_integ_ary[1] - logl_integ_ary[0])*(theta_integ_ary[1] - theta_integ_ary[0])*(logM_integ_ary[1] - logM_integ_ary[0])\n",
" integ = 0\n",
" for logl in logl_integ_ary:\n",
" for theta in theta_integ_ary:\n",
" for logM in logM_integ_ary:\n",
" integ += integrand_Ncumul([logM, logl, theta], J_fid)\n",
" integ *= measure\n",
"\n",
" return integ\n"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "f86656e00a8f4f01a8513acbd9aa1e83",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HBox(children=(IntProgress(value=0, max=30), HTML(value='')))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Ncumul_manual_ary = np.array([Ncumul_manual(J) for J in tqdm_notebook(J_ary)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plot results"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x2ad72fd2b4a8>"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAhQAAAFfCAYAAAAWDVXXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl4U1X6B/DvSbojEALINqKk4gKuofx0hnFGJVVBKgJdUBRxkFTUcVwbcRlxXLB1R0UbVBRE6OIy1hmXVmccZRyHNuiouDYoM1I22xQQKF3O74/clCQkaZqlN8v38zx5NPee3Lz3kPS+OdsVUkoQERERhUOjdgBEREQU/5hQEBERUdiYUBAREVHYmFAQERFR2JhQKIQQBiGEWe04iIiI4hETioMsAHRqB0FERBSPmFAAEEIYATSoHQcREVG8YkLhpAfQrHYQRERE8SohEgohhE4IYRZCVPnZbxZC5CuPEq99JillXd9ESkRElJhS1A4gXEp3hQHOFgaDj/1mAM1SymrluUEIUS6lLBZC6MCWCSIiorCJRFl6W0kslkspJ3htb/CxrVFKma20VtiVzblwdn2Us8WCiIiod+K+hSIQpQXC6GOXQ+nqKHMrqwegYzJBRETUewkxhiIAAwCHj+0e3SNK60YugFwhhKmPYiMiIkoYCd1CAf+zNxxwW3NCSmkDUNCbA2dmZkqN5mA+NmTIEAwdOhStra0YOHBg9/bePN+xYweGDh3amzAC8n6vcMsH2u9rXzDb3J9Hsy78xRNO+UjXRzx/NgKVCXY7vyv8rvjblkifjUBlIvFd8d4Xbn3s2LEDO3fu7H6+d+/e/VLKTJ+FpZQJ8YCza6PBa5sJQKOPslUASsJ5vyFDhkhfFixYEPLzCRMm+DxmqLzfK9zygfb72hfMNvfn0awLf/GEUz7S9RHPn41AZYLdzu9K4G3J/F1JpM9GoDKR+K5474t0fQDYIf1cFxO9ywNwtlJ4C3tFzNTUVJjNZtTU1Hhsz8vLC+t5JPX22D2VD7Tf175gtrk/j2ZdhHL8vq6PeP5sBCoT7HZ+VwJvS+bvSiJ9NgKVicR3JVp1UVNTA7PZDADt/sok9CwPZVBmi5RSeJVtAGCRYQzAzMnJkfX19SHH6+eYiPQx4xXrwhPrwxPr4yDWhSfWh6dI14cyczLH176EbqGQUjoA2JXEwl1MzuZQsj8C68Ib68MT6+Mg1oUn1oenvqyPRGqhMAEolYeuOWGGM4EoU54bARRLKYvDeb+xY8fKs846C3l5eVFvfiQiIlJTTU0NampqsHz58u+klGN9lYn7hEIIYQCQD+e0TxOAMjgHYlrdypjhXMBKB8Ag3dafCFU0ujyIiIhiWaAuj7ifNiqltMOZRPhNEtyTCyIiIoq8hB5DEU2tra0+Z3kQERElGrdZHn4X2Yj7Lg+1sMuDKL7s2rUL27dvR3u731lvREknJSUFGRkZGDp0KDIyMnosn9BdHkREPdm1axe2bduGUaNGITMzE0KInl9ElOCklOjo6MCePXuwefNmDBs2rFergnpjQkFECW/79u0YNWoUsrKy1A6FKGYIIZCamopBgwYhPT0dW7duDSuh4BiKEHEMBVH8aG9vR2am79sPEBGQmZmJtrY2v/s5hiKKOIaCKH58+eWXOP7449UOgyimBfM9SdqVMomIiKhvMKEgIiKisDGhCBHHUBARUbIIZgwFE4oQDRw4EFarlffxIKKYYbFYMGjQIEyYMKHHckIIWCwW2O32Q/aXlZXBYrGgrKwMVqsV1dXVAIDq6mo4HI6oxE6xLS8vD1arFQBa/ZXhtFEiogRRWlqK7Ozs7kTBYDD4LWswGFBaWuqxzW63o7i4GMXFxSgpKene7nA4YLVaUVpaioaGhoAxWCwWOBwOlJeX9yp2h8MBnc77xtDq8BdLqOeWLNhCQUSUQPR6PQoLC7tbFbzZbDZMnDjR576CggJYLBbk5+d7bNfpdCgsLPTZmuGtqKgIxcW9v5lzZWVlr18TLf5iCfXckgUTCiKiBFNcXOz3V3Rzc7PPX99lZWXQ6/UwmUw+X6fT6Vx96AEZjUYYjcbeBQygtra216+JFn+xhHpuyYIJRYg4KJOIYpXRaERzczNsNlvQrykvL0dBQUHAMrm5ueGG5pOrKyEWxFIssSSYQZkcQxEi16BMIqJYZDabUVFR4fGLuq6uDiaTCXV1dYeUt9vtyMnxuV5RN++uEF/HcHUJ1NbWwmazYcGCBTAYDFi0aBGam5vhcDiwfv367vEbroGedrsdZWVl3bG7WlHKyspgNBrhcDjQ3Nzs0UpitVqh1+tht9u7z7O8vByLFi3CggULkJOTg4KCAtjtdtTW1qKqqqr7PV3xGo3G7laZQLF4n5t3DAA84gvm3ONJXl4e8vLysHz5cg7KJCLytmBZs9oh+LT8Kn3YxygqKsLkyZODungFMzYiGK6BnhaLBYCzpcT1XKfTdV/0KyoqupOb/Px8GAwG2O12j4GggHNMx6JFi7pfZ7FYUF1djfz8fBQXF6OgoAAmkwkOhwMTJkxAY2Njd3Lgel+LxYKcnJzuVoeysjKP98nNzYXBYIDBYAgYi/e5ueIrLS3tHvzqcDiQm5uL2traoM490bDLg4goARmNRuj1eo9uD9cvaW+BZoP0lvf4DL1eD4fD4fEerot2IHa7HTabzaOFpaioqHtsSGVlZXeLius9vc/V9b46na47QVi/fr3HgFWj0eizxaanc7PZbIfMpNHpdDAYDN2t16Gee7xiCwURJa1ItATEMtfgzPLy8h5/FRsMBtTX1wccdNjTVFR/QpkOWldXB51O53Gxd3VHuOJ1H2Da3Nx8SGy+YnV1ewDO83E4HCGNmaivr/d5/OzsbI+ptbEyFbYvsIWCiChB5efnBz0ds7i42ONi60tvBnm689cy4o/rQm8wGGAymbof+fn5aGxs7I63vLy8e42MRYsWHXLx9nUxdzgcKC4u7m5F6OmC7681IVAS0tx8sCutt+cez9hCESO6ZBc0gvkdEUWOa2xAdXV1jy0LJSUlqKioOKSbQQ2uGCoqKvyWca23UV9fD5PJFHTLyZgxY7Bp0ya/yYavWHwd22Qy+YyvsbExarNhYh2vYCGK5LTR1v27cObTF+Al2+sRiIyIktn69es9nhcVFWHJkiVBJQlVVVWwWCw+xxRYrdYeZ3n44/6L3cX94u09rkCn08FkMh0yBgRAd2yumRaBkgnvBMF1LPdkwlXGbrd3d5t4x+KL0WiEwWDwiM/hcKC+vt5jJkpP5x4vgpk2ql28eHGfBZRIVq5cubimpgbHHnts2Me6960X8c+mV/G3TW/jh52tOPvoX0Gr0UYgSiICgJ07d2Lo0KFqhxF1rvtvNDU14YwzzkBGRgbGjx+PtrY2TJo0CYBzauRjjz2GhoYG7Nu3DxMnTkRGRgYAYNCgQZg7dy6qq6uxbt06rFu3Dl988QUaGhqCWtTKbrfDYrHgww8/hF6vhxACS5YsQUNDAzQaDSZNmoTq6mosXboUjY2NGDx4MMaNG4eMjAy0tbVh3bp1aGpqwrRp0wAAc+fOxapVq/DFF19g06ZN2LhxY/e+zMxM5OXlYeXKlSgrK8PSpUuxZcsW5ObmwmazwWKxoL6+HhqNBuPHj0dGRgZGjBiBlpYWbNiwAfv27YPD4UBhYSEqKyuRlZWFSZMm+Y3F+9wmTJiAwsJCrFq1Ck1NTdi4cSPq6uqwdOlSZGRkwGazBXXusSTQ9+TYY49FXl4e7rrrrk2LFy/2uWaCkFJGNcBElZOTI+vr6yNyrL1tXbi26ln8fesDkOjEWF0OVs1+AkMPS/w/gER94csvv8Txxx+vdhgUITabDXV1dYdM7SwrK0NjYyPvtRGiYL4nQogGKaXPBUvY5REDstI1sM65AtecvAJpYgi+ddTD9Gwe/v3f0AZAERElsvLycp/dLyUlJQk7JTMeMKGIERohcF3uJDx+3ssYpD0Fe9p3YE7FRVj+8SqwFYmI6KDc3Fyf4zzq6uqSdkBkLGCXR4gi2eXhbdP2fVhYfS++3bsGAJBruBCP5t2LjNSMqLwfUaJjl0fiqaur81gXw263Q6/XhzxwlMLv8uC00Rg05vBMVF3+J9xYfSL+tu0u1Npfw9QVX+OFwqdwhO4ItcMjIlJdIi5dHe/Y5RGj+mdqsOziQlx14mpkil/gh11fYsrz0/G3xvfVDo2IiOgQTChiWIpW4PpzT8WDpmoM0f4a+zpaccWr8/HwB0+gS3apHR4REVE3JhRx4LyTh2JVkRXHZTpvnfvkx49gXuWV2N22W+XIiIiInJhQhCiSK2UG45iR6Vh12U3IPXwpUtAf6/77LqY8Nx1f7/i6T96fiIiSVzArZXKWR4iiOcsjkAMdEkvf/gqrvroee+S3SNNk4oGp92PacdP6PBaieMFZHkQ948JWSSYtReDGqcfhvrPWYLh2Cg507cMf3vgD7n73PnR0dagdHhERJSkmFHFICIHzjYNgnfUwxmfeDAEtnt/wLC5eMw/New+9EQ0REVG0MaGIY+NHp+GZSxYgd2g50qBHQ9NHmLLiAny+7XO1QyMioiTDhCLOHT5Qi0cvPhNXHL8WAzQnYue+JsxaXYiXP39F7dCIiCiJMKFIAOmpAjdMNeCeM1/AqJSZ6OhqQ8lbN+O2t+/Egc4DaodHRH3EYrFg0KBBmDBhQo/lhBCwWCwxdTMtm82G3Nxc3o8jTjGhSBDOcRUD8fSs+3By1u0QSMXaz15E4UuXYMfPO9QOj4j6QGlpKUpLS2G323tMFAwGA0pLS7vvhRELjEYjSktLgyprsVhQXFwc0vs4HI6QXhcNvmIJ59zUlPQJhRDCJIQwCiHyhRAlascTrnFHpMJ6yVxMHfYM0sXh+GxbA6asyINtC2+FTpQM9Ho9CgsLUV1d7XO/zWbDxIkT+ziq4Ol0uqDKFRUVhXzRraysDOl10eArlnDOTU1JnVAIIXQAyqWUNillNYBiIUTspOshGjJAiwcvmoTicWuh05yKlv07MHvNRVj9yUu8FTpREiguLkZ5ebnPfc3NzUFftGOZ0WiE0WgM6bW1tbURjiZ0vmIJ59zUlNQJhZTSIaXM9toWOx2KYUhLEbj2vNG47+wVOCJlNjplB/5YdwdufnMR2jra1A6PiKLIaDSiubkZNhtbJr1ZLJaY6fKIpVgigbcvVwghzACC67yLE0IInHtKfxw1dDFurxmHT/bci1c3VuHL7V9h+cxlGDlgpNohElGUmM1mVFRUePzSraurg8lkQl1dnc/XuLpJ7HY7jEZj9y3CbTYbFixYAIPBgEWLFqG5uRkOhwPr16/vHvPgKqPX67t/dVutVpSWlsJisbiWbQ74PsGy2+3dXQK1tbVBxed6X4fDAbvdjrKysu56crXYlJWVwWg0wuFwoLm52SNmq9UKvV7fHTMAlJeXY9GiRViwYAFycnJQUFAAu92O2tpaVFVVBTxXf7E0Nzd7nJs7VwwAPOIL9vyjLSESCqXrohBArpSywMd+MwDXik8GKWWZ137XpzkhV4U6dlQqyufMxt2vj8U7W2/EVzs/w/nPT8fTM57AaUecpnZ4RKrJfjC750IqaLypMexjFBUVYfLkyUFfUMrKylBScnAYWW5uLgwGAwwGQ/dgSYvFAp1O131Braio6E5SjEYjli9fDovF0n0Ms9mMhoaGoN8nWK4Bpa73CiY+AMjPz4fBYIDdbveIAQAKCgqwaNGi7tdaLBZUV1cjPz8fxcXFKCgogMlkgsPhwIQJE9DY2Nh9XNd7WywW5OTkdLc6BDpXf7HodDqPc3OPz30QrcPhQG5uLmpra4M+/2iL+y4PIYQRgAnOZOCQT6QrmZBSVivjJKqFEB6di1LKOimlFcAiIUR+X8Td1/T9tSgrmoirT1yDQZrTsOtAMy6pnIsXbKs4roIoARmNRuj1eo9uD9evW1/Wr1/vMZDTaDR6tGTo9Xo4HA6PC7/rghiI93iNnt4nWN7HDTU+wNl6YLPZPFpzioqKusehVFZWIicnx+N9vevV9d46na47QQj1XL3PzWazwW63e5ybTqeDwWCA1WoN+/wjJe5bKKSUNgA2JbHwpVhKOcGtvN3VIqEkGxOklK7htHYAEwH4Hh4d51JTBBbmjsJxI5bjnvcewA/tq/Cn9xbj861f4t5zFyNNm6Z2iER9KhItAbHMNTizvLy8x1+qriZ6wHmBdTgch/TvR2IwZzDvE6pQ46urq4NOp/O42Lu6IwDnhdl9MGtzc/MhLSq+Wlgida719fU+j5+dne3RAqT2YNu4b6EIROkK8ZVoOJSkwg7AvbXCAKCiL2JTixACZ5/YD9bC25Bz2D3QIB2vbKxA/otzsPPnnWqHR0QRlJ+fH/QUSYfDgeLi4u5fvL4uToFaOIIVzPuEKpT4XBd6g8EAk8nU/cjPz0djozPhdCVmDocDVqsVixYtOiRuX+fR23P115oQKAlpbj7YUx+Jf59wxH0LRQ8MAHz9SzTDOZbCqqw/YVDKlistHj3asWNHdxMY4OwrdB/AE+uOHpGK8kuLcPefx+DNpuvxxQ4bpqyYjmfzn8JJw09SOzwiigBXf311dXWPYxTGjBmDTZs2+b0wBkOn03lc4Hy9NhLvE0muro6KCv+/JV1re9TX18NkMgU93qO352qz2Xwe22Qy+YyvsbEx6quKWq3W7oRIMcRf2YRuoQCgh++Blg4AOgBwja2QUpYp4yiCMnToUNTX13c/4imZcNH106C06DRce3IFBmpORvP+rShYPRuvfP6a2qERUYjWr1/v8byoqAhLliwJuK6BazyA+4XPdcGz2+3dSYJ3suBeDjjYj+/O1QLQm/cJVU/xAYeOK9DpdDCZTIeMNwHQ3QXimm0RKJnwfp9gztVXLL4YjUYYDAaP+BwOxyHXnmDOv7fMZrPHtQ6A36Zs7eLFi8N6s1hx1113jQBwweLFi61u2wwA8hYvXrzUq2whgC2LFy9eF+r7Wa3WxfGYRHjTaARyDANx7ICpaPh+K1q7vkDtd+/AsXcvzjjqV9CIRM85KRns3LkTQ4cOVTuMqLNYLCgrK0NTUxPOOOMMZGRkYPz48Whra8OkSZMAOKcrPvbYY2hoaMC+ffswceJEjBkzBi0tLdiwYQP27dsHh8OBwsJCVFZWIisrC5mZmViyZAkaGhqg0WgwadIkVFdXY+nSpWhsbMTgwYMxbtw4ZGRkoK2tDU1NTdiyZQtsNhv0ej3Ky8sxatQoTJ48OeD7jBgxAhaLBR9++GH3+/hit9u7y+n1egghgooPQHeM69atQ1NTE6ZNmwYAmDt3LlatWoUvvvgCmzZtwsaNG7v3ZWZmIi8vDytXrkRZWRmWLl2KLVu2IDc3FzabDRaLBfX19dBoNBg/fjwyMjIwYsSIgOc6adIkn7F4n5vrviyFhYVYtWoVmpqasHHjRtTV1WHp0qXIyMiAzWYL+vwDCeZ7ctdddzW5X2fdiUQZ4a8MylzuPgBTGSdRJaUc5FW2FkCt9/TR3hg7dqw866yzkJeXh7y8vJDjjiX2re24+bXn8emeByDRiZwRZ8A66zEMzBiodmhEYfnyyy9x/PHHqx0GxSGbzYa6urpDppmWlZWhsbHR74qk8SjQ96SmpgY1NTVYvnz5d1LKsb7KJPrPz3ooXRte9ADCWkJu4MCBsFqtCZNMAIBheCpWzJ2PvFFPIxU61Dd9gCnPzcC3O79VOzQiIlWUl5cjP//Q1QRKSkpi6k6t0ZaXl+caS9Hqr0xCJxRSSgcAuzLbw51OStn7ic9uWltbYTabUVNTE85hYs6ALA0eKDwL1526FoeJsdi29wdcsHIW3vnmXbVDIyLqc7m5uT7Xjqirq0uq26zX1NS4xmv4bbJOpC4PE4BS9y4PZbsZzgSiTHluhHNtirBu5ZaTkyOVASoJ6+8bW3DbO7dga0cdAIGFE6/Djb+5GkIItUMj6hV2eVA46urqPBaWstvt0Ov1Plsu4lkw3xMhRIOUMsfXvrifNqpM+cwHkAvAKIQoBdDomrGhTA01KwmHDs7povF3X1gVnDluEF4a+gSufflxfL5nGZ5a/wg+3/YVll1Yhqy0LLXDIyLqE321dHW8i/suDymlXZnymSulFFJKi/f0TymlVVleuzqcgZjuErXLw9uRQ1Oxet51OH/Uo9CiHz7Y/CamPJeP/7X+T+3QiIiojyRVl0dfS4YuD3ddXRLlf/sCT37ye+yTm5GVMgjlM57Er47kzcUo9rHLg6hn4XZ5xH0LBfUNjUZg4eQT8PiUCgzWno69HS2YWzUXz/77RbVDIwoKfzwR+ReJ7wcTCuqVs8YdjoqLn8UxWZdCogP3/eNO/OH1W3Gg84DaoRH5lZKSgo6ODrXDIIpZ7e3t0Gq1YR2DCUWIkmUMhS9jhmWgYt4fcebQu6FBGt74pgIzVl7Cm4tRzMrIyMCePXvUDoMoZu3atQv9+/f3u59jKKIo2cZQ+NLRKfHAmx9j5dfX4YDcAV36CLxQUI4Tho9XOzQiD/v378fmzZvxi1/8ApmZmZz6TARnN0d7ezt27dqFlpYWjB49Gunp6QFfE2gMBROKEDGhcJJS4tWG/+Gef/werV2fIUVk4P5zSzHjhGlqh0bkobW1FT/99BPa2trUDoUoZmi1WvTv3x96vb7HZAJgQhEVTCg8ffrDHvzhz7fjvwecXUCXnbIQt0++gTcXIyJKIJzlEQXJPIbCl5OPPAwVcx/EhIE3AtDghU+ewpw1Zuxu2612aEREFCaOoYgitlD4tv+AxB2v1eH1zSXowC6M6JeNVUXlGKMfo3ZoREQUJrZQUJ/JSBMoLTDh5pw16CfGoOnnRkx7YQb+Zv+H2qEREVEUMaGgiNMIgSvOPA6PT63A4Sm/wf7O3Vjwynw8+dGzXFyIiChBMaGgqPnt8YOx5uJyHJc1HxJdeHjdfbj6tZvR1sFR9kREiYYJBUXVUYenYc28W2A6vAwaZODtxldxwQuzsW3PNrVDIyKiCGJCESLO8gjegCwNnpwzEwuOX4l0MRzftfwH5z07HbYfN6gdGhERBYGzPKKIszx6T0qJ121N+NP718LRtQFakYo/me7F7JNnqR0aEREFgbM8KCYIITB9wkg8n78SR6bno1O247baEtz61t3o6OKNm4iI4hkTCupzJ47OQsXc+/B/ulshoEXF58+jaPXlaN3fqnZoREQUIiYUpIqhA7VYMfd3KDhqOVIxCJ9s+yfOe/ZCfLPzG7VDIyKiEDChINVkpAncO+s3KJm4FoeJY7B932ZMXzkLb39Tq3ZoRETUS0woQsRZHpGhEQK/++3RWHb+WgxPycWBrr246vWFeOgfT3ARLCKiGMFZHlHEWR6R98P2dlz78uP4/OdlACR+e+R5eGL6A8hKy1I7NCIiAmd5UJw48vBUrLrsOpw7/BFo0Q/v//AWpj2fjx9bf1Q7NCIi6gETCoopA7I0WHrRNFw5fjUyxRH4YdfXmPL8hfh487/VDo2IiAJgQkExJ0UrcP15J+D+syug156On9ubcUnlpVhlW6N2aERE5AcTCopJQghMO3UYVsx6BmPSL0YXOrD4vdth+eudXASLiCgGMaGgmHbC6Ey8NHcxTtfdCYFUVG98EUWr58Gxz6F2aERE5IYJBcW8wwdqsfzSS5A/ejlSoccn2z7ClBUX4tud36odGhERKZhQUFzIStfgvvxf48YJa3GYOBbb9/4X01fOQt1376odGhERgQkFxRGNRmDBWdl45LyXMCzFhLaun1H8WjEeX/cUF8EiIlIZE4oQcaVM9Zw9XocXCh/HcZkLAUg8+tGDWPja9djfvl/t0IiIEhJXyowirpSpvp92d6KkugYf/HQ7OrEPY/Un4vmCpzG8/3C1QyMiSkhcKZMS0uD+Wjx5yXRccvRKZIiR+Lb5M0xdcSE+2fKJ2qERESUdJhQU1zJSBW6ffipuPX0tdBojWg/sQOGai/Dy56+qHRoRUVJhQkFxTyME5kwahSfOfw6/SJ2JTnkAJW/dhLvfLUWX7FI7PCKipMCEghLGL4/tj5UXLcEph90CAS2e32DFvMor8fOBn9UOjYgo4TGhoIRy5OEpePbS+Thv2BNIQX+s+++7yHuhEFt2bVE7NCKihMaEghKOrp8Gj1yUi/nHvYhMMRo/tH6FqSsuRMOPNrVDIyJKWEmfUAghTEIIsxCiVAhhVjseiozUFIGbzx+Pu89Yg0Gaidjd/hMuWjsHL3/2Z7VDIyJKSEmdUAghDAAgpbRKKS0ASoUQRpXDoggRQmDG/w3H8hnPYnTaLOdgzbdvwL3vPcTBmkREEZbUCQUAI4Bit+d1AHwu2EHx69Qx/bB6zn0wDigBoMFztmW4vPJq7D2wV+3QiIgSRkgrZQohTgFgAjAYgA6AHkAzAAeAnwDUSSnjYnUhIYROSulQ/r8BgEVKWdfT67hSZvz5eX8Xbnv1Hfz1xxJ04mccOWAcXrpoOVfWJCIKUkRWyhRCnC2EqBRCfAtgEYAhcCYRdQCsyn+ble23CiHWCyGeEkKcHfYZ9BybThkHUeVnv1kIka88Stz3uSUTBgDNwSQTFJ/6ZWjwUNG5uObk1cgQo/DDro0477kL8cmW/6gdGhFR3OuxhUIIMQZAKYBGAJVSyg1BH1yIUwEUARgD5y//70MP1e97GAEYlKeLpJQTvPab4UwUqpXnBiWWYq9y5d7bAmELRXx785NtuOO9a9DSZYNWpOP+c8ow88RpaodFRBTTArVQBEwohBCT4ezauF9K2RpGAAPhbNV4R0r5XqjH6eE9jACW+0goGnxsa5RSZrs9LwFgdbVWBIMJRfzb+L+fcfWrd2Bzm3Pmx/xTr8Wis6+FEELlyIiIYlNIXR5KywSklIvCSSaUY7RKKW9xHlYcFc6xekMIoYNz4KU3hxDCpJQxAah26/ow9VV8pK5xv+iHqsvKkKO7AYDAsxuW4vLKa3kbdCKiEPhNKKSUm6SU70byzaSU70aj2yMAA5wDRb01AzAorRpVAGqFEC1CCKns69GOHTuQk5PT/bBarZGLmvrMkAEpWHnZVZh55GPQIgsf/PevOH/FbGzfs13t0IiIVGe1Wj2udXCOk/QpmDEUZ0c295G7AAAgAElEQVSrmyKSfHV5KK0N5e7dG8r2KgDrpZRlob4fuzwSS5eUeOb9z/BYw1XYL5swIHU4VhU9gxOGH692aEREMSPcWR6WCMdDFHM0QsB85klYel4VdNqTsKt9K2atLsAbX9aqHRoRUVwIJqHIFUK8JYSY35fjHyJI72ObLtyDtra2wmw2o6amJtxDUQyZPH4E1sx+EUekT0GH3Ic//GUhHvnHcoSyXgsRUaKoqamB2WwGgIH+ygTT5eG+RrGEc0xCLZzrTtT18ZgIv/x0eegAtEgphVfZoBew8oddHontp92dKF77KDa0LgMAnGMowNLpdyNVm6pyZERE6gm3y6NOSqmBc0nqRQBsAArhXMyqUQjxrRBimRBiphBiQMSijgBl5oZdSSzc6biAFQUyuL8Wq+fdgGlHlEKDdLxjr8KFL1yG1v1hTXgiIkpYwSQUDgCQUtqklGVSylwlwTgHwIMAWgFcCaAaQIsQYlnUog3MV9cG4FyUq/suokpLRtjJBLs8El96qsAjhbNwzckrkAY9vmr+GLnPzEDjT5vUDo2IqE9FqstjjJQy4F9QZeEqE5xJxgIA9VLK/+t1xCFQVr7MB5CrxFAGoFFKaXUrYwZgh3PshCGc2R0u7PJILjUbNuGOvy3E7q5vka4ZgKemL8Nvs3+pdlhERH0q5JUyQ3wzHYBKOJfpfiaiB48hTCiSz2ebHbjy1euxtf0fENDiljP+hCtOm612WEREfSYiNwcLlpTSIaU8Bwl+G3B2eSSfE0fr8Mpl5Rh32KWQ6MSSD25DyV/uRWdXp9qhERFFVUS6PEIlhLhZSvlAVA4eA9hCkbz2HZD4Q+VKvLf1Xkh04tRhZ+GFosfQL62f2qEREUVVn7VQCCHeFkLcpzwdHMljE8WKzDSBpy+ei9+NK0cKBmDDtr/h3GcK8GPrj2qHRkSkmkh3eQwBcIsQYgmca1UkLHZ5JDeNRuDWqWfh7t+uRZYYjaa9X2PKilnY8ONnaodGRBRxqnR5CCFOBWAP9w6lsY5dHuSy3r4DV79+NX7qaECKyMBDUx/FtONz1Q6LiCji+npQ5oZETyaI3E00DEXlnOdxZHoeOuR+/OEvC/HEume5XDcRJZWIJxREyeiooVmouuxBnDzgKgASj3x0H27+y53o6OpQOzQioj7BhIIoQgYPSMGquddj8uFLoEEaXv1qNeasWYA9B/aoHRoRUdQxoQgRB2WSL/0yNHjy4gIUjSlHKnSob/oHLni+EFt2bVE7NCKikKm6DkWi46BMCqSrS2JZ3dco/+wq7JU/YGD6UKwsfAYnDDtB7dCIiELWp4Myicg5rfTq3GPxx0kvQaeZgNa2HchfPRu13yb0bGoiSmJMKIiiRAiBgtOH45HznsWIlPPR3rUPV/55IZ5dv4IzQIgo4TChIIqy34zrj6dmPICx6VcCkLjv/XtwR+1izgAhooTChCJEHJRJvXHikel4qvB6TDjsHgikYs1/XsT8ajNngBBRXOizQZlCiLellOeGfaA4wkGZFIrtrZ247eUP8EHzDWhHK8bqj8OKgmcwov8ItUMjIupRVAdlCiFmAZgohJgf7rGIEt3hA7V4aPZvcMHIF5AlRuPb5q8wY1U+vtn5jdqhERGFJRJdHrcAmADgyggciyjhDcjS4E8FJ2BO9ioM1JyCHXu3omD1bDT82KB2aEREIQsroRBCTAbQIKXcBKBBCHF2ZMIiSmwZqQI35Y3ClSeVY4j2t9jT3oo5FZei7rs6tUMjIgpJuC0UJQDuV/6/FIAlzOMRJY0UrcCCyYNx66THMDJlBtq72nDlawtR8WmF2qEREfVayAmFEGIMnIM6vwcApZWiVQhxSoRiI0p4QghM/7/+KJtyDwypV0CiC7fW3orH1j3JtSqIKK6E00JRgkNbJJYAWBTGMeMGp41SJP3y2Aw8kX8TTsi6BYDA0o8exq1vL0ZnV6faoRERRW/aqBBiIIB3fU0dEULUAzhbSrmr1weOI5w2StHQ1NKJkpdfw0eO2yDRjjOPOg/LLnwY6SnpaodGRBSVaaO3wNka4UvStFIQRdqIQVo8efEMTBm+DFr0w9+/fwuzX7ocu9t2qx0aEVFAoSYUuVLKl33tULbnhh4SUXIbkKXBQ7PPwpzs55CGwfjP9o9xwQtF2L5nu9qhERH51euEQghxBYCehqGXCyFuCi0kIkpLEbjjwgm47tTVyBKjsXnX1zh/RT42NW9SOzQiIp9CaaEollI+EKiAlHI5gOLQQiIiANAIgeLJY3HvWasxQDMezW0/4oKVBdjw43/UDo2I6BC9SiiUhazeDbL4y0KImb0PiYjcXWAciacveAFDtKdjb0cLZq+9GLXfvK92WEREHnqVUEgp35VS3hJk2VuklK+EFhYRuTvt6EF46eJncET6VHTIfVj4+gK8aHtN7bCIiLrx9uVEcSJ7WCaqL3sY4w6bC4lO3PnejXjsg+fUDouICAATipBxYStSw5ABqai4/A78csj1AIClH9+LP779EFfVJKKoCmthK2Vp7Z8iuUCVEGIAAL1rue54xoWtSE3tHRLXVr6I2i13Q6ITU7IvwmPT74JWo1U7NCJKYCEtbKXcm6MoUncQFUKcCsCcCMkEkdpSUwSemH0JZo15CBqk4c3GNbis4loc6DygdmhElKQCdnko0z+zhRBPhXrTLyHEUUKIpwFMllI+GMoxiOhQWo3AkpnT8LtxT0OLfvjox7eQv2o+fj7ws9qhEVESCupeHsq9O0oBTABQB6AWQDMAu3uXiNKlYQCgB3AOABOA9QBukVK2Rjx6FbHLg2KFlBLlf9+ARxuK0Y5mHDngRLx86QoMyhykdmhElGACdXn0+uZgStdFEZzJggHOARoOADrlv3YA9QDq/C3PnQiYUFCsWfuvr3H3uiuwX27B4ZkGVF/yPEYNHKV2WESUQCKaUJATEwqKRW/957+w1M3Hnq5GDEgbjsqLX8DYIUerHRYRJYho3G2UiGLQeScdgaemrYZOcxJ2HdiKGatmw8aluomoD/SYUCT6Tb6EEDohhFkIUap2LESR8KtjhmJF/vMYmvIr7OtswcVr5+DvjevUDouIElwwLRRFUY9CXa6mG52qURBF0EmjB+Kli6z4Rdp5aJd7seDV+Xjt87+qHRYRJbBgEooJQogKIcRRUY5FFVLKOjgHkhIlFMOwTKy99FEcnVmELrTjprf+gOfWr1E7LCJKUMEkFA4AVjgTi/uFEEuEEFckaoJBlEhGDErFS3Pvxkn9r4REF+59/3Y89P5TaodFRAkoJYgylVJK1y3LXwa616UwCSGuBCABNMI5TfT7qETZAyGEDkAhgFwpZYGP/WY4180AAIOUsqwv4yNS0+D+WqyceyMWrhmEj5rvx7L1D6Ktcz8WnXUdhBBqh0dECaLHhEJKeaWPba1wJhfeCYYFzkWt3gHwbl8kGEIII5zrYTQr//XebwbQLKWsVp4bhBDlUsriaMdGFCv6Z2rwzCXzcc3afvjb9j/iWdsTaOtow+JcC5MKIoqISE0bHQPnQlfFAAoALFf+P+qklDYlWfA3DqLYlUwo5e1wxkqUVDLSBJ6YPRuTD78fAlq8+J/luPWtu9Alu9QOjYgSQFgJhTKWYj2ABjgTCBucF3CNlHJRJAIMh9IVYvSxyyGEYFJBSScjTeCx2TORO+whCKSi8otVuPkvtzGpIKKw9TqhEEKcotwsrBPKYE04WySypZQ5yg3FYoUBzkGl3rq7R5TEohhAjtI9QpTQMtMEHi6ahnOGPQIN0vHaV5W49vWb0NHVoXZoRBTHghmUCcDZGgHnhdcIQMDZGlEeYwmENz0ODsZ057r3iGvaaF1vD7xjxw7k5BxcfdRsNsNsZj5C8SEzTeDBwnNhqUzD29v+gDe//TPaXzuAJ6Y/glRtqtrhEVGMsFqtsFqt7puG+Cvb4708hBBPATDDmUQAzlaJUinlpjDjjChlcOZyKeUEt20mOJOebK+yVQDWhzPbg/fyoESwt60LlsoP8fa2a9CJn/Gboybj6QsfR3pKutqhEVEMCvdeHq6xEWZlbMSVsZZM9EDvY1vYq2K2trbCbDajpqYm3EMRqSYrXYP7C36NqcPLkYIB+Mf372L+y8XY175P7dCIKIbU1NS4WuEH+isTTAtFM5ytEgDwHVRcbyIQPy0UOgAtUkrhVbYBgEXp7ggJWygokfy8vwu3v9yAN5sWoh0tmDDyNKzIX45+af3UDo2IYki4LRR1UspbpJS3AKhCkCtmKmMuVCWldACwK4mFO104yQRRoumXocE9syZg6ojlSBND0LDlY1xaOQ+723arHRoRxYlgEoolrv+RUrZKKV9WEoxFCJxgWCIfbkC+ujYAoBTOMSAAulsywk4m2OVBiaZfhgZ/mnUypg57BuliOD7dasPFay+FY5+viVJElEwi0uXRG8qKmYVwJhNjpJTaiB3c/3saAOQDyIVzwaoyAI1SSqtbGTOcC1/pEKGlt9nlQYlq974u3PXy13hzWzH2yx8xdvBxeLHwBQzp53dwNxEliUBdHhFNKLzetFlK6a/VIO4xoaBEtntfF+5+pRFvbi3GXvkDxgzKxuqiVRh22DC1QyMiFYU7hiJUCT1GgV0elMj6Z2pw+4xsTBn+DPqJbGxqaUTRmouwZdcWtUMjIhX0eZdHMmELBSWDXXu7cM8rP+DtbVdjj/waI/v/Ai8VvYgjdEeoHRoRqUCtFgoiinMDsjS4beaRmDKiHAM047Fl9/9QtPYibGqJp6VoiKgvMKEgooAGZmlw24xRmDbqKQzUnIxte5pQtOZifPfTd2qHRkQxhAlFiDiGgpJJ/0wNFl04AjNHL4NOk4Of9m5H0ZqL8PWOr9UOjYj6AMdQRBHHUFAyamuXePQv21FtvwHNXf9C/zQdXix6AScMO0Ht0IioD3AMBRFFRHqqwA3TDsecYx7DYM2vsfuAA3MqLsWnTZ+qHRoRqYwJBRH1SmqKwDVTBuPycY9iqPZM7DmwC3MqLkX9j2yxI0pmTChCxDEUlMxStALF5+iw4KQHcbj2HOzr+BmXVV6Of23+l9qhEVEUcAxFFHEMBRHQJSVe+mA3nmq4HVs7/4JUTTqWzyzHGUedoXZoRBQFHENBRFGhEQJzzuiP60+7DyO009He1YYrXlmA9xrfUzs0IupjTCiIKCxCCOT/6jBYzrgbo1IK0NHVjitfW4i3v31b7dCIqA8xoSCiiMjL6Yc7zrwTR6RcjE7ZgWte/z3e+OoNtcMioj7ChCJEHJRJdKjcUzLxp9xbcWTKPHTJTlz3xvV49YtX1Q6LiMLEQZlRxEGZRP796+v9uOOtx2Bvt0JAg4emPojp46arHRYRhYmDMomoT51+bAbK8q6HIc0MiS7c+Neb8MaX7P4gSmRMKIgoKk41pOGhC66HIfUKSHTh+r/cgL98/Ve1wyKiKGFCQURRc9JRaXho+g0Yk3Y5uuAcU/HWN5z9QZSImFAQUVSddFQ6Hpl+M8akzUWX7MDvX78W73z7jtphEVGEMaEgoqg78ch0PHbhIhyVOgdd6MDVr/8edd+9q3ZYRBRBTChCxGmjRL0zfnQaHp95O45MuwhdsgNX/flqvPfd39QOi4iCwGmjUcRpo0Sh+ep/B7DwlcXYfKACWpGG5TOexm8Nv1U7LCIKAqeNElHMOO4XaXh61mKMTs1HpzwA86sL8Y9NH6odFhGFiQkFEfW5Y0el4emCP+GItJnokG1Y8IoZH37/T7XDIqIwMKEgIlUcOzId1vx7cETqdHTINsx/eQH++cO/1A6LiELEhIKIVHPMyHQsL7wfv0jLQ4fcj9+9fAX+9cO/1Q6LiELAhIKIVDV2RBqeLSjFqNSpaO/ah3kvz8fHmzngmSjeMKEgItUdPSIdK4oexMi089DetRfzqn+H9f+zqR0WEfUCEwoiignZw9PxQtFDGJl6Dg50/YzLKi/HZ1s/VzssIgoSEwoiihmGYRlYOfsRDE89C21de3Dx2nn4Zue3aodFREFgQhEirpRJFB1jhmVg+azHMET7S+ztaEHh6kvxfcv3aodFlNS4UmYUcaVMouhqsLei+LUr0NJlgz59JP58WQVGDhipdlhESY0rZRJR3JlgGIhHp1oxQHMCmtu2IP/FS7Hz551qh0VEfjChIKKY9evjBmHJZCsOE2Oxbe/3yH/xUjj2OdQOi4h8YEJBRDHtvJOH4Y5fP4MscST+u/sbFL00D7vbdqsdFhF5YUJBRDEv/7Rf4IacZ5EhRuK7ls8wZ+0V2Htgr9phEZEbJhREFBfm/eYoFJ/wDNLF4fhiRz3mVV2Jto42tcMiIgUTCiKKC0IIXJN7DC4duxyp0KOhaR2KX/k92jvb1Q6NiMCEAkIIkxAi3/VQOx4i8k+jEbj5/PEoOOpppGAAPtj8Lq6tuQmdXZ1qh0aU9JI+oQBQLKWsllJWAygSQujUDoiI/EvRCtw+3YhpI5+EFll457s3UPLmbeiSXWqHRpTUkjqhUFokmt02rQdgUikcIgpSeqrA3bNOh2noY9AgHa99WYXFdXeDC/URqSepEwoABgDuk9odACaqFAsR9UJWugb3F5yJM/QPQiAFqz9diQf+8bDaYRElrYRIKIQQOiGEWQhR5We/2W2cRInX7p+8nrPLgyhODMjS4IGCc3HagPshoEX5+mV48qNlaodFlJTiPqEQQhjh7KZohrPFwXu/GUCz2ziJaiFEuVuRwV4v4TJ8RHFkcH8tHsy/AKf2uwuAwMPrHsJz9SvUDoso6cR9QiGltCmJgt1PkWJlv6u8HQfHSXi/RgfnOAoiiiMj9Fo8MKsAJ2beBgC49+/34KVP16gcFVFyifuEIhBlxobRxy6HEMKkJBrurRrZAOr6JDgiiqijDk9B2YxLcHzGTQCAP9begVe+eEXlqIiSR4raAUSZ96BLF/fukXJltocDQK2UMqgujx07diAn5+AdXM1ms+te8USkkmNGpqJs+hW44dU2fHvgcZS8aUFGSgamHjtV7dCI4pLVaoXVanXfNMRfWZEo06yUsRTLpZQT3LaZAJRLKbO9ylYBWC+lLAv1/XJycmR9fX3I8RJR9Hyy6QBuev1hbGpfDo1IQfmFT+Hs7LPVDoso7gkhGqSUOb72JXSXBxElp1PGpOHeKddhdMql6JIdWPjaVfjw+w/VDosooSVDQqH3sS3sqaGtra0wm82oqakJ91BEFAWnHZOBP51zC0alFKBDtmPBK8X493//rXZYRHGppqbG1a0/0F+ZRO/y0AFokVIKr7INACxSypAHYLLLgyg+1H66F3e9dxuaOl9HhjYLq2evwikjTlE7LKK4lLRdHsoAS7uP+3PowkkmALZQEMWL3JOzcMsZd2OY9lzs79yLuZWXY+P2jWqHRRRXkq2FwgSg1L2FQtluhjOBKFOeG+Fcm6I4nPdjCwVRfKn+Zyse/PgG7Oj8OwakD0LlRWswdshYtcMiiisJ3UIhhDAoy2lbABiFEKVKEgEAkFJaoaw7oUwPNYWbTBBR/Jn1ywG4xvggBmt+hV1tLbho7aXY1LJJ7bCIEkbCtFD0NbZQEMWfLinx3HvNeOrTq+DoqsfQrBF4eU4FRg0cpXZoRHEhoVso1MIxFETxRyMEfneWHr8b9zgGak7Cjr1NmL12Drbu3qp2aEQxLanGUPQ1tlAQxa+OTomlbzZh5Tdm7O76EkcONKDy4jUY0s/vIoBEBLZQEBF5SNEK/H7KCBQalqGfOBo/tNpxaeXl2N22W+3QiOIWE4oQscuDKL6lagVuPP8IzBz9NDLFEfjmp4244pUr0dbRpnZoRDGHXR5RxC4PosSw74DEndVfombLPBzATzBln4Nl05+AVqNVOzSimMMuDyIiPzLTBG6/8DicPeRJpOAw1DW+g9vf+SP4Y4uod5hQEFHSG5ClwV0zjfjVwEegQToqP1+LRz58VO2wiOIKEwoiIgBDBmhxz8zfwNjvPgho8eTHT+CFhhfUDosobjChCBEHZRIlnlGDU3DvhedjfMatAIC7/3Y3ar7id5yIgzKjiIMyiRLXZz8cwE1/fhLfHXgCWpGCZ2c9gzOOOkPtsIhUx0GZRES9cOKRaVh87kIckXIxOmUHil9diE+bPlU7LKKYxoSCiMiHXx6bgdvOWoTh2qlo69yHy6p+h8afGtUOiyhmMaEgIvIj9+Qs3HD6PRismYTdBxyYUzEPTbub1A6LKCYxoQgRB2USJYeZpw9A8ckPYYDmROzYuwVz1s6DY59D7bCI+hQHZUYRB2USJY+uLonH3vwRK76ah5/lJpxw+KlYM3slstKy1A6NqE9xUCYRURg0GoFrzhuF/KOeQroYjs+3b8CVr12D9s52tUMjihlMKIiIgpCqFSjJOxrnj3gSqRiIdZvfx01/taBLdqkdGlFMYEJBRBSkjDSBP844CWcNWQotMvHG13/G3e/ey/t+EIEJBRFRr/TP1OCeWb/ELwc+AIEUrPzkeTzx0TK1wyJSHRMKIqJeGtxfi/tnnQNjv7sBCDz6z4exasOLaodFpComFCHitFGi5DZCr8X9M2difMYtAIC73l2M1ze+rnJURNHBaaNRxGmjRAQAn28+gBtfewLfHXgSGpGC5TPKcabhTLXDIooKThslIoqSE0an4a7zrsLolEvQJTuw8LWrUf8jf2xQ8mFCQUQUptOPycBtZ92CEdoLcKBrP35XvQBf7fhK7bCI+hQTCiKiCDCdnIkbfnUXhmrPws/tu3BJxWX4vuV7tcMi6jNMKIiIImTG//XHlac8gEGaiWjZvxNzKuZi255taodF1CeYUBARRYgQAnN/q8Pc4x5Ff804bN3zIy6puIw3E6OkwISCiCiCNBqBhecOQ+GYJ9FPjIG95VvMq56PvQf2qh0aUVQxoSAiirBUrcAN5x+B80c+hQwxAp9t+wTFry1EW0eb2qERRQ0TihBxYSsiCiQjTeDWC7NhGroMqdDjn5s/xA1/uRGdXZ1qh0bUa1zYKoq4sBURBeOn3Z24pXI93m8xoxM/o/DEItx3zr0QQqgdGlGvcWErIiKVDO6vxeIZOfi/wx6FBumo/KwCD3zwgNphEUUcEwoioigbpU/B4uln4JTM+yGgRfm/y/Fs/bNqh0UUUUwoiIj6wNEjUvHHaedhfPqdAIAlf1+Ct755S+WoiCKHCQURUR858cg0LDonH9mp10BC4ro3boBti03tsIgiggkFEVEfOv3YdFz362KMTJmJ9q42XPGyGT84flA7LKKwMaEgIupjU4yZmH/y7dBrfonWthbMq5qPln0taodFFJakTyiEEDohhFkIUap2LESUHIQQuPS3A1CU/SAOE2OxuXUTFrxyJRe+oriW9AkFANd8Wp2qURBRUtFqBK6ZOgxTRj6OdHE4NjTV46a/lqBLdqkdGlFIkj6hkFLWAbCrHQcRJZ+MVIFbLsjGGYOWQot++Os3b+ChDx5SOyyikCR9QkFEpCZdPw3uvNCICf1KIaDF0/9+Gmv/s1btsIh6LaYTCrfxDVV+9puFEPnKo6Sv44s0q9Wqdggxg3XhifXhKdHqY6ReizsvMOG49FsAAHfU/hHvb3o/qNcmWl2Ei/XhqS/rI2YTCiGEEYAJQDMAg4/9ZgDNUspqKWU1gGohRLn7fiFEiZ+Hqc9OpBf4RTiIdeGJ9eEpEevjuFGpuP3cOTgy5XJ0yU5c9edr8OX2L3t8XSLWRThYH56YUACQUtqURMHf+IZiZb+rvB3OBMT13CqlLPPzqAs3vtbWVp/bve8+2tvnkdTbY/dUPtB+X/uC2eb+PNp3bo31+ojnz0agMsFu53cFOP2YdFz/6+sxTHsu9nfsxbzq+Wja3XTI65Ltu5JIn41AZSLxXemDO2D7vdtozCYUgQghdACMPnY5+qr1weFw+NyeyF8EJhTB72dC0fvt/K44TZuQhStOvhcDNadi595tmFc1H7vbdif1dyWRPhuBysRJQuF3RmTM375c6fpYLqWc4LXtXSnlIK+ytQCqpJRBt/EoCUgxnN0q5cG+VgixD4D7/K4dAHbCmb25N1/05vkQ5RiR4v1e4ZYPtN/XvmC2uT+PZl34iyec8pGuj3j+bAQqE+x2flcCb0vm70oifTYClYnEd8V7X7j1MQTAULfnGillpq+CKWG8iZr0cI6t8OZAL9eTULo/et0F4q9CiYiIklFcdnkQERFRbInnhELvYxtXuyQiIlJBvCYU9fCdPOgB8F7AREREfSwuEwoppQOAXZnt4U4XiSmhRERE1DvxkFD46toAgFIAZtcTZeYHkwkiIiIVxOy0USGEAUA+gFw4F6wqA9DoPq1TWS3TDmf3h0FKWaZGrJGktLoUAsiWUlq89pXA2aVj6M3U2Hjmrz6EEOVSymL1IlNHgPowwjn1WQ9lBVmVQuwzAeoiHwd/iNiTpdVSmQJvAJANt7+Vyvbu1txk+GwAAevD79/YRNXDZ+OQ7aGK2WmjysqXZcrDX5lEvKj6vJ26sqx4qZTSLoQoF0IYlDpKdIfUh5JsmoUQhW7l6qWUuX0amTp8fj4A5Lj9kTAnyefD32ejSEpZoDyvQhK0XCrnDbfPQIsQol5KaYNzVeHu+hBC1Cndxgmrh/rw9x1KSP7qAs5lFvzVUUjiocsjqfi6nbrygei+QEgpi5PgYgHA7+3lDQAGSSkHKYubFSiPhOfn86GD5/nr4HudloTi57Nh8trWrLRYJDojnAv0udQByFHO3f2zsB5utyhIYD7rA/D7uUlk/urCbx2FKmZbKMiDEQeXFdcB0Cdo60xQ3JuwXQNzE/0XVyBSSocQwi6EaARggbOZP2nrA56/PB3wcXPBRCOlrBZCuLfEGOC8aBqh/BJVOABMBJDQ3R4B6iPp+KsLKWVdpOuICUUUuPXR5bqaGr32m3HwV0MwYz9cLRR1yuurwm2a6ktRqA935ngbOxOl+rAAWK484qZfOAp1UczLzgcAAAnDSURBVAfnuCsXI+JoKnk49eFKIpUWzWblgmEE8JPXYeKmqT/S9dEHIUdNNOoi0nXEhCLC3AbH9XjbdeW5IYgBhg44195wcd1ZNeb/UEapPlyvjbtfntGoD+UPjVlKWaDUSa0QIuYHI0ajLpQxRhVKU79deTRG5QQiLIL1YfEaTzTYa39ctF5FsT7iTh/URUTqiAlFhCmtBjblA+BLsfuNzpQ/gD31adrh+asrbkSpPrpfC2efcNyIUn0UQmnCVsrnwlk3MZ1QROuz4T6LQQihh2cyHrMiUR/KTDD3Fio7PBMKHeLkOxOl+ohL0ayLSNYRB2X2IRHibdeVX5ruzZQGxPjFIhih1ocbIxKoXzSM+rDj0F8tcXHR8CfUuhBC6JSZHd3ipWswkGDqQ/lvtVsztklJrtw/G9lIkr8dvuqjD0PsM+HURaTriAlF3zLAd3NjdzOWOHg79RylGcvFIoQoUbbVJsIfSYRXHy5x0XwbpJDqQ0k4DUKIfGWbMQHWGgi1LhwAKoQQJiFEia++5jgVsD6UX65VcHZ3tQghJA72p5crnw0TnH87EuE7E3J9BPE3Jd6EVBc9fGZCwi6PvtXjbdeln9upu5q8ohpd3wu5PpR9cdkNFEA4n49Em/UTTl24kqm4/yXuJmB9KH8fBvl6YayPpQlRuPWRSHUScl0E2B4StlAQERFR2JhQ9D3edt0T68MT6+Mg1oUn1ocn1sdBMVEXTCj6Fm+77on14Yn1cRDrwhPrwxPr46CYqQsmFH1I8rbrHlgfnlgfB7EuPLE+PLE+DoqlumBCET287bon1ocn1sdBrAtPrA9PrI+DYrouYvb25fFKJOlt1/1hfXhifRzEuvDE+vDE+jgoXuqCCQURERGFjV0eREREFDYmFERERBQ2JhREREQUNiYUREREFDYmFERERBQ2JhREREQUNiYUREREFDYmFERERBQ2JhREREQUNiYUREREFLYUtQMgovii3NWwEM57BkwEUJ5sd3hMNMp9ILIBVMB5A6oCKWWxulFRvGFCQUS9tUhKaQG6k4sWIUS2lNKuclwUOh2cN58qgfMulUwmqNfY5UFEQVPueqhzPZdSOgBUgxegeOeQUmZLKYWUMpfJIYWCCQWRyoQQtUKIBuWh6/kVqiv0EadBlUgoqoQQpcrnsiVOPpukIiYURBEghDAIIaqU5EAqj1ohhCmY10spJygPh49jlyjHqhJClCsPk7JPp/R/9ybWEuUiIYUQjUKIkh7KSuWCUiqltEspB3nFaQSwvjcxqM3t34sXSYUQIl95mF2fLymlRUo5Ac5uEKKAOIaCKAKUJuICABBCSAB1UsrccI4phDACqAJQ7X0sIYRJCJEPIBdAYy9jLQNQJoRoUY5dFqisECIbQKmvZnAlRgQ6hlqU7hmDnwGjOgAmOAcgHpLERVsPsamhHoDdlSgqiabPBJfIH7ZQEEWQ6wILZyIQ7nEaABS7BkC6Uy5EdgC9ap3wUhnk6xsD9KmXApgQRgzR5LcbRkppU1pa1BorEFNdREp9uCcPNgCL1IqH4hNbKIgiy9XFEe4vzyoAZYF+wUopbUIIaxjvUQ7ALIQw+Xsfpenb375SOBOeWP0VW4AwE7soimpsSjdYT4lerZSyWilvlFLa3PY1I8aSHop9TCiIIisXzhHzIf/yVS4GBgBLgiheBecYhl5TEhIbnDM0/CUuRl/dGUqMFa7z9HFBUpWSCJkRgwlFX8QmpQw60VRaw94FMMhrV3NEg6KEx4SCKLJMcE6jDEcB3PqzA5FS1vkaWKhcJExwdovoAej8jHMoB1AuhNB5v59yXF+DRF3HtStl9ABy4Gwm90kptwie4z0q3d8zUMzKvuVwJloFyktc750LYIFb/79rbAkAWIQQrvKlUkpXzO+6jqXUob/jG+BcvMuivNaMg7/eB/vqjurhPALGFkJdTFb+61qMqtfjdpTEcoHXZpNbnETBkVLywQcfEXjA+UdYAjD38nW1Xs9bvLeFEIf3MUvgXNHSu6xOibnExz4znBcy920Gpbz3w+85K+/RCOcgRI+YehuzK1b3YylxHlJfynuaAsTV4r0/0PG9z9HPtmDPw29sIdSFTnm0eP979eIzY1COZYZzXIx3vVSFemw+kuehegB88JEoD+UPsfS+cAbxOu+LhwRQFUYcPi9WynEPuSgoF4tGH9sPSTJCjKcKzl/g7ttMzt8zvYtZKed9YTW6H6unY7rtb/Bx4fR7fB9lS73/nXp5Hv4Sit4co9TXMSL9YELBRzAPzvIgihwTnF0V4c4ccMBtNcreUJrDDXBOA/Rmh7Nrwls5AIPbDBXXcSI1JiIfzl/z7uqhrK7Zy5gdcCYC8NoWKYGO7x3fT+5PQqx7DyEcI67W/6DExjEURBGg9MsbAYQz68KlHgFG2CtjGAqUMq5+9jrpvJmT63UmIYT3Sy3wcaGSzjEEdjjHOLj69E0yAmtLuCUpHkmWdI53cNVVb2OO9mBBn8eXPY9p6XXdR+AYsTrDhpIQEwqiyHBNF/X+JR6KKjgHShp8tXZI5xTPOqB7Ea1qeXBwoEMp09uBoeVwNuGrIdSYe8XXwNMIC/k83GLr7TE4E4NiBrs8iHpJWbbaO3FwjYgPe+VD6ZzyZ4fzF2kw3Jve64HulRh7w6q8zqzMRIjIxV0enEoaKJ5QY+6twigfP5zzcMXWV3VBFHFMKIh6QWnC9xhfoHR3mAFYI/gLuADKolO9eZE8ePfPfO99Qgij+zgJH69z3bbaZ8tIGKpxsCvFPZ7/b+/+bhqGgTiO/7xBVTaADWAE6AYVGyBGQIwQNmCGsgEZATpCNmAG9+EuUkjd4MSJCuL7eapaxTqnD774z+XS61dMijnT6sTn2U3ox1FsC98LYFEkFMA4jayCZXf2YCdp73sYZuFP9jeypY/qRK2JSuk9G0+SHhNPuXdxuPhUJUuW5i5H/SDbE9AfDLedeHJjbmtP5KhlNSRa/WSvPW7Z/y63/YvE9bn9GIott431iFiBxYUY47ljAP4UnzXYyJYariR9xhGVCRPtvceBgkT+NtCNbNBpZMcFV/IZkVSVyk4hqS+/Zp0T40+xTOXxVB57Ix3vExiK2QfXSvbk3siOdr748syzLBGqZcco606br/JTGwNtVe21ue17u/ey/+FN3wtrZd37VGyZ9+Laf9vKZso+5kxmU0IIO3X6CKSQUABnttQgDsyFhAI5WPIAAADFSCgAAEAxEgoAAFCMhAIAABSjUiZwfvsQQvv+iFs2vuG38KPJo2qh4P/ilAcAACjGkgcAAChGQgEAAIqRUAAAgGIkFAAAoBgJBQAAKEZCAQAAih0ACnbvx9I4UwcAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 576x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(8,5))\n",
"\n",
"plt.plot(J_ary, Ncumul_ary[:,0], label=\"MC integration\")\n",
"plt.plot(J_ary, Ncumul_manual_ary, label=\"Manual integration\")\n",
"\n",
"plt.xscale(\"log\")\n",
"plt.yscale(\"log\")\n",
"plt.ylim(1e-1,1e4)\n",
"\n",
"plt.xlabel(\"$J$ [GeV$^2$\\,centimeter$^{-5}$]\")\n",
"plt.ylabel(\"$N(>J)$\")\n",
"\n",
"plt.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 2. Differential SCD calculation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The same as the cumulative case, but during integration $\\delta N$ is phase-space-constrained to lie between $J$ and $J + \\delta J$:\n",
"\n",
"$$\\frac{dN}{dJ}(J) = \\frac{1}{\\delta J}\\int_{J < J(M,l,\\theta) < J + \\delta J} dM\\,dl\\,d\\theta\\,\\rho_M(M)\\cdot 4\\pi r^2 \\rho_V(r),$$\n",
"making sure that the final result is not sensitive to the specific choice of $\\delta J$. You may have to increase the number of points by quite a bit to get a smooth SCD!"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"# def integrand_dN(x, J_fid, dJ):\n",
"# \"\"\" Integrand for calculating differential SCD by integrating \n",
"# over the mass function, distance from GC and opening angle from GC and to the Galactic plane \n",
"# \"\"\"\n",
" \n",
" \n",
"# M200, r, theta = np.exp(x[0])*M_s, np.exp(x[1])*kpc, x[2]\n",
" \n",
"# # If J-factor is within than specified values, add the corresponding \"dN\"\n",
"# # (impose phase-space constrain)\n",
"# l = np.sqrt(l**2 + R_sun**2 + 2*r*R_sun*np.cos(theta))\n",
"# if J_fid < J(M200, c200_SC(M200), l)/(GeV**2*Centimeter**-5) < J_fid + dJ:\n",
"# return 2*np.pi*np.sin(theta)*pref*r*M200*rho_M_r(M200, r)\n",
"# else:\n",
"# return 0."
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {},
"outputs": [],
"source": [
"# npoints = 5e6\n",
"\n",
"# def dNdJ_MC(J_fid):\n",
"# \"\"\" Perform cumulative SCD integral with scikit-monaco\n",
"# \"\"\"\n",
"# # Set value of differential element\n",
"# # Check to make sure the calculation is not sensitive to this choice!!\n",
"# dJ = 1e-1*J_fid\n",
"\n",
"# dN_val, dN_err = mcmiser(integrand_dN, npoints=npoints, xl=[np.log(M_min/M_s),np.log(R_min/kpc),0], xu=[np.log(M_max/M_s), np.log(R_max/kpc), np.pi], nprocs=5, args=[J_fid, dJ])\n",
"# return dN_val/dJ, dN_err/dJ"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [],
"source": [
"# dNdJ_ary = np.array([dNdJ_MC(J) for J in tqdm_notebook(J_ary)])"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {},
"outputs": [],
"source": [
"# deg = 180 / np.pi\n",
"# sr = 4 * np.pi\n",
"# srdeg2 = sr * deg**2"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [],
"source": [
"# plt.plot(J_ary, J_ary**2*(1/srdeg2)*dNdJ_ary[:,0]/(GeV**2*Centimeter**-5))\n",
"# plt.xscale(\"log\")\n",
"# plt.yscale(\"log\")\n",
"# plt.xlabel(\"$J$ [GeV$^2$\\,centimeter$^{-5}$]\")\n",
"# plt.ylabel(\"$J^2\\,dN/dJ$ [GeV$^2$\\,centimeter$^{-5}$\\,deg$^{-2}$]\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment