triNFWtutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# importing\n",
    "import numpy as np\n",
    "from math import *\n",
    "import matplotlib.pyplot as pyplt\n",
    "from matplotlib.colors import LogNorm\n",
    "# cosmology stuff.\n",
    "# If you dont have cosmolopy installed, run:\n",
    "# git install cosmolopy \n",
    "import cosmolopy.distance as cosdist\n",
    "# triaxial NFW\n",
    "from triNFW import triNFW"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up the cosmology in cosmolopy\n",
    "cosmo       = {'omega_M_0' : 0.27, 'omega_lambda_0' : 0.73, 'h' : 1.0}\n",
    "cosmo       = cosdist.set_omega_k_0(cosmo)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('halo.psi:', 1.0471975511965976)\n"
     ]
    }
   ],
   "source": [
    "# construct a triNFW object\n",
    "mass      = 1E14   # Msun/h, since h = 1\n",
    "concen    = 3.0    # concentration\n",
    "zd        = 0.4    # cluster redshift\n",
    "overden   = 500.0  # 500 times of `wrt`. If overden < 0.0, then it is caluclated as virial density\n",
    "wrt       = \"crit\" # overdensity type: crit or mean\n",
    "qa        = 0.5    # the axis ratio, 0 < aq := a/c <= qb := b/c <= 1.\n",
    "qb        = 0.8    # the axis ratio, 0 < aq := a/c <= qb := b/c <= 1.\n",
    "theta     = pi/7.  # First Euler angle\n",
    "phi       = pi/9.  # Secon Euler angle\n",
    "psi       = pi/3. # Third Euler angle\n",
    "\n",
    "# construct this halo\n",
    "halo  = triNFW.pNFW(mass    = mass,\n",
    "                    concen  = concen,\n",
    "                    zd      = zd,\n",
    "                    overden = overden,\n",
    "                    wrt     = wrt,\n",
    "                    qa      = qa,\n",
    "                    qb      = qb,\n",
    "                    theta   = theta,\n",
    "                    phi     = phi,\n",
    "                    psi     = psi,\n",
    "                    cosmo   = cosmo)\n",
    "# check psi\n",
    "print(\"halo.psi:\", halo.psi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "triNFW/triNFW.py:216: RuntimeWarning: invalid value encountered in sqrt\n",
      "  case_smaller_than_one   =   1.0 / (1.0 - X**2) * ( -1.0 + 2.0/np.sqrt(1.0-X**2) * np.arctanh( np.sqrt((1.0 - X)/(1.0 + X)) ) )\n",
      "triNFW/triNFW.py:217: RuntimeWarning: invalid value encountered in sqrt\n",
      "  case_larger_than_one    =   1.0 / (X**2 - 1.0) * (  1.0 - 2.0/np.sqrt(X**2-1.0) * np.arctan(  np.sqrt((X - 1.0)/(1.0 + X)) ) )\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Sigma')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEOCAYAAACaQSCZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3Xd8leXdx/HPL4uQAGHPsAPIkBnZMtQqKEvEAYKjyCqgrfZ5Wvto97JWKwqIIIgTpDjAiaLMAEoYsiHsBBTC3iPkev5IaFPKOISc3Gd8368XrzZ3zjn51t7my3Vf931d5pxDRETkSiK8DiAiIsFBhSEiIj5RYYiIiE9UGCIi4hMVhoiI+ESFISIiPlFhiIiIT1QYIiLiExWGiIj4RIUhIiI+ifI6QEEqW7asq1GjhtcxRESCyrJly/Y558pd6XUhVRg1atQgNTXV6xgiIkHFzHb48jpdkhIREZ+ERGGYWXczG3/48GGvo4iIhKyQKAzn3EfOucEJCQleRxERCVkhURgiIuJ/KgwREfGJCkNERHyiwsj1yarv2XXopNcxREQCVkg9h5Ff+46dZvg7ywGoV6E4na4rR+d65WlRvRTRkepUEREAc855naHAJCcnu/w8uOecY0vmMeZsyGTOxr0s3X6As+ccxWOj6FC3HLfUL0+nuuUpFR/jh9QiIt4ys2XOueQrvk6F8d+OnjpLyuZ9fL1hL3M2ZpJ59DQRBsnVS3NLg/L8qEFFapaNL4DEIiLeU2EUkOxsx6pdh/lq/R6+XLeHDT8cBaBO+WL8qEEFbmtYkcaJCZhZgf5cEZHCosLwk/QDJ5idWx7fbDvAuWxHpYRYbmtYkdsaVqRlzdJERqg8RCR4qDAKwcHjZ/hqw15mrf2B+ZsyOZ2VTZn4GG5tWIHbr69E61plNGkuIgFPhVHIjp/OYu7GTD5b8z1zNuzl+JlzlIqL5raGFbmjcSXa1CpDlMpDRAJQ0BWGmdUC/g9IcM71yT3WCfgDsBaY6pybe7nP8LIw8jp19hzzNmXy6ervmb1uD8fPnKN0fAxdG1Wke5PKtKxRmghdthKRAOFrYfj1OQwzmwR0A/Y65xrlOd4FGAVEAq865/7qnNsKDDSz6Xk+wgHHgFggw59ZC1JsdOS/5jROnT3H3I2ZfLL6e95fvou3v9lJxRKxdG9SiZ5Nq9CwcglNmItIUPDrCMPMOpDzC/+N84VhZpHAJuBH5JTAUqCvc25d7ven5xlhRDjnss2sAvC8c+7+y/28QBlhXMqJM1nMXr+XmSt3M2/TXs6ec9QuF0/PplXo1bQK1crEeR1RRMJQQIwwnHPzzazGBYdbAptzRxSY2VSgJ7DuIu/Pzv2vB4Ei/ktaOOJioujRpDI9mlTm0IkzfLr6Bz5cuYvnv9zE819u4oYapbizWSJ3NK5EQtFor+OKiPwHL5YGqQKk5/k6A2hlZmWAPwHNzOxJ59xfzKw3cBtQEhh9sQ8zs8HAYIBq1ar5NXhBKhkXQ79W1ejXqhq7Dp3kwxW7+GDFLn71wWp++9FaftSgAn2aJ3JjnbKaLBeRgOD3Se/cEcbHeS5J9QG6OOceyf16ANDKOTfiWn9WoF+SuhLnHGt2HeG95RnMWLmLgyfOUr54EXo3T+Tu5ERqlyvmdUQRCUEBcUnqEnYBVfN8nZh7LOyZGdcnJnB9YgK/ur0+X2/Yy/Rl6UxYsJVx87aQXL0U99xQlW6NKxEXo3UjRaRweTHCiCJn0vtmcopiKdDPObf2Gn5Gd6B7UlLSoLS0tGvOHGj2Hj3FB8t38W5qOlszjxMfE0mPplXo17Ia1ydqW1oRuTYB8RyGmU0BOgFlgT3Ab5xzE83sduAFcm6rneSc+1NB/LxgvyR1Jc45UnccZOq36XyyejenzmbTqEoJ+rWsTs+mlYkvolGHiFy9gCiMwhbqhZHX4ZNnmbFyF+98s5MNPxylWJEoejWrzP2tqlO/Ugmv44lIEAmrwgj1S1KX45xjRfoh3l6yk49X7eZ0VjbJ1UsxoE11ujSqSJGoSK8jikiAC6vCOC+cRhgXc+jEGaYvy+CtJTvYvv8EZYvF0K9lNe5vXZ0KJWK9jiciAUqFEcaysx0LNu/jjUXb+XrjXiLN6NKoIg+3q0nzaiW1FImI/IewKoxwviR1JTv2H+fNxTt4NzWdo6eyaJKYwMPtanJH40pael1EgDArjPM0wri046ezeH95Bq8t2s7WzONUKFGEB9vWoF/LapSM017lIuFMhSEXlZ3tmJeWycQF21i4eR9FoyO5JzmRge1rafFDkTClwpAr2vDDESYu2MaHK3dxLtvRpVFFhnSoTZOqJb2OJiKFKKwKQ3MY12bPkVNMXrSdt5fs4MipLFrXKs2QjrXpVLecJshFwkBYFcZ5GmFcm2Ons5j67U5eXbCNH46con6lEgzrVJvbG1XUirkiIUyFIfl2JiubGSt3MW7eFrZkHqda6TiGdapN7+ZV9CCgSAhSYcg1y852fLFuD2PnbmZVxmEqlohlcIda9G1ZjaIxKg6RUBFWhaE5DP9yzrEgbR+j52zm220HKFsshsEdanF/q+pa8FAkBIRVYZynEYb/fbvtAC99ncaCtH2UiotmUIdaPNCmBsVUHCJBS4UhfrV850Fe/CqNuRsz/1UcD7apoRGHSBBSYUihWLHzIKNyi6N0fAxDO9ZiQOsamuMQCSIqDClUK3Ye5B+z05i/KZOyxYowvHNt+rasRmy0ikMk0KkwxBNLtx/g+S82sXjrfionxDLy5jr0aZGohQ5FAlhYFYbukgo8KZv38eysjaxMP0SNMnE8fms9ul1fiYgIPTkuEmjCqjDO0wgjsDjn+Gr9Xv7+xUY2/HCUBpVK8L9d6tFRS46IBBRfC0PXCcRvzIxbGlTgk0dv5IV7m3LsdBYPvbaUfhO+4bv0Q17HE5GrpMIQv4uMMHo1q8Lsxzvyux4N2bTnKD3HpDD87eVs33fc63gi4iNdkpJCd+x0FhPmb2XCgq2cycqmf+vqjLwpiTLFingdTSQsaQ5DAt7eI6d44as03l2aTlx0JD/pnMTD7WroVlyRQqY5DAl45UvE8uc7r2fWT2+kZc3SPPP5Bm5+bh4zVu4ilP4iIxIqQqIwzKy7mY0/fPiw11EkH5LKF2fiQzfwziOtKBkXzWNTV3Ln2EUs23HQ62gikocuSUlAOZfteG95Bs/O2kjm0dN0b1KZX3a9jioli3odTSRk6ZKUBKXICOOe5KrM/XknRt6UxBdrf+Dm5+byjy83cfLMOa/jiYQ1FYYEpPgiUTxxaz2+eqIjt9SvwKiv0rjpubl89N1uzW+IeESFIQEtsVQco/s1Z9qQNpSOj2HklBXcO34J63Yf8TqaSNhRYUhQaFmzNDNHtOfPd15P2p6jdHtpAb+esYbDJ856HU0kbKgwJGhERhj9WlVj7s87M6B1dd5asoPOz83l3aU7yc7WZSoRf1NhSNBJiIvmdz0b8fHIG6ldLp5fvLea3i8vYs0u3VYt4k8qDAlaDSqXYNqQNjx/TxMyDp6g++iFPP2hLlOJ+IsKQ4KamdG7eSJfPdGJB9vU4O1vdnDz83N5f3mG7qYSKWAhURh60lsSikbz2x4NmTmiPYml4nh82nf0nbCEzXuPeh1NJGToSW8JOdnZjqlL0/nrZ+s5efYcQzvWZnjnJC1qKHIJetJbwlZE7t1UX/+8E90aV+alrzdz2wvzWZCW6XU0kaCmwpCQVbZYEf5xb1PefqQVEWYMmPgtP3t3JfuPnfY6mkhQUmFIyGuXVJbPHruRkTcl8fGq3dz8/Dz+mZquSXGRq6TCkLAQGx3JE7fW49NHbySpXDH+Z/oq+k/8hh37tUWsiK9UGBJW6lQozrQhbfhjr0asSj/MbS/M55V5W8g6l+11NJGAp8KQsBMRYfRvXZ0vH+9Ihzrl+MtnG7hz7CLW7tZt2SKXo8KQsFUxIZZXBrRg7P3N+f7wSXqMTuHZWRs4dVb7bohcjApDwpqZcfv1lZj9eEd6Na3CmDlb6PbSQpbv1PawIhdSYYgAJeNieO6eJrz+45acOJ3FXS8v4g8fr9MufyJ5qDBE8uhYtxxfPN6R+1tVY+LCbXQdNZ9vtx3wOpZIQFBhiFygWJEo/tjret4Z1IpzznHv+MX8duZaTpzJ8jqaiKdUGCKX0LZ2WT5/rAMPtK7O5EXb6TpqgUYbEtYCpjDMrJaZTTSz6RccjzezVDPr5lU2CV/xRaL4Xc9GTBnUmuzc0cbvP1qnO6kkLPm1MMxskpntNbM1FxzvYmYbzWyzmf0SwDm31Tk38CIf8wtgmj9zilxJm9pl+PyxDvRvVZ1JKdu4fdQC3UklYcffI4zJQJe8B8wsEhgDdAUaAH3NrMHF3mxmPwLWAXv9G1PkyuKLRPGHXo14+5FWnDp7jj4vL+KZzzdwOkujDQkPfi0M59x84MKLvi2BzbkjijPAVKDnJT6iE9Aa6AcMMrOAuYQm4atdUlk+/1kH+rRI5OW5W+g5OoX13x/xOpaI33nxC7gKkJ7n6wygipmVMbNxQDMzexLAOfd/zrmfAu8AE5xz/7Xgj5kNzp3jSM3M1H4HUjhKxEbztz5NePWBZPYdO0OP0QsZM2cz57K1Aq6EriivA5znnNsPDL3E9yZf5n3jgfGQs+OeX8KJXMItDSrwZfVSPPXhGp6dtZGvN+zl+XuaUL1MvNfRRAqcFyOMXUDVPF8n5h7LN+3pLV4qFR/D6H7NGHVfUzbtOUrXUQt455ud2m9DQo4XhbEUqGNmNc0sBrgPmHktH+ic+8g5NzghIaFAAopcLTOjZ9MqzPppB5pVK8mvPljNI6+nknlUu/tJ6PD3bbVTgMVAPTPLMLOBzrksYAQwC1gPTHPOrfVnDpHCUrlkUd78cSt+3a0BCzbvo8sL85m9bo/XsUQKhIXCsNnMugPdk5KSBqWlpXkdRwSATXuO8tOpK1n3/RH6tqzG093qExcTMNOGIv9iZsucc8lXel1I3KaqS1ISiOpWKM4Hw9sypGMtpi7dSbcXF7Iq45DXsUTyLSQKQyRQFYmK5Mmu9Xn7kVacPHuO3mMX6fZbCVoqDJFCcH4hwy6NKvLsrI30nbCEXYdOeh1L5KqERGHotloJBglx0bzUtxnP3d2EtbsO0/WF+Xyy6nuvY4n4LCQKQ3MYEizMjLtaJPLpYzdSq1wxhr+znP+d/h3HT2uvDQl8IVEYIsGmepl4/jm0DSM6J/HPZRl0e2khqzM0QpbAFhKFoUtSEoyiIyP4+W31mDKoNafOnqP3yylMmL+VbE2IS4AKicLQJSkJZq1rleGzx27kpuvK86dP1/PQ5KV6QlwCUkgUhkiwKxkXw7j+Lfhjr0Z8s3U/XUctYEGaVl+WwKLCEAkQZkb/1tWZOaI9peKiGTDxW/762QbOnvuvVf1FPKHCEAkw9SoWZ+aI9vRtWY1x87Zw97jFpB844XUsEd/XkjKzRuRsqRp7/phz7g0/5boqWktKQtUnq77nl++tAoNn7mrM7ddX8jqShCBf15LyqTDM7DfkbJfaAPiUnP24Fzrn+lxjzgKVnJzsUlNTvY4hUqDSD5xgxJQVfJd+iPtbVePpbg2IjY70OpaEkIJefLAPcDPwg3PuYaAJoFuSRApB1dJxTB/ahiEdavH2NzvpNSaFLZnHvI4lYcjXwjiZu592lpmVAPbyn7vmiYgfRUdG8OTt9XntoRvYc+QU3V9ayAcrMryOJWHG18JINbOSwARgGbCcnI2RRKQQdb6uPJ8+diONKifws3e/4xfTV3HyzDmvY0mYuOoNlMysBlDCObfKH4GuheYwJFxkncvmhdlpjJm7mbrlizPm/mYklS/udSwJUgW+gZKZNTazHkBzIMnMel9LwIKkpUEk3ETlLivy+sMt2XfsND1Gp+gSlfidr3dJTQIaA2uB808ROefcj/2Y7apphCHhaM+RU4x8ZwXfbj/AfTdU5bc9GuouKrkqvo4wfN1guLVzrsE1ZhIRP6hQIpZ3BrXi+S83MXbuFlamH2Ls/c2pVa6Y19EkxPh6SWqxmakwRAJUVGQE/9vlOl57OOcuqh6jU7Q5kxQ4XwvjDXJKY6OZrTKz1WYWcJPeIuGuc73yfPLojdStkLM5029mrOF0lu6ikoLh6yWpicAAYDX/nsMQkQBUuWRR3h3Shmc+28CrC7exMv0QY+5vTmKpOK+jSZDzdYSR6Zyb6Zzb5pzbcf6PX5OJSL5FR0bwVLcGjOvfnK2Zx7njxYV8vWGP17EkyPlaGCvM7B0z62tmvc//8Wuyq6DbakUurkujSnw0sj1VShblx5NT+fusjZzTjn6ST77eVvvaRQ7rtlqRIHHq7Dl+M2Mt76am07Z2GV7s24yyxYp4HUsCRIGuVhssVBgilzctNZ2nP1xDybhoxvRrTnKN0l5HkgBQoM9hmNmLFzl8GEh1zs242nAi4o17kqvSqHICw95exn3jl/Cr2+vzcLsamJnX0SQI+DqHEQs0BdJy/zQGEoGBZvaCn7KJiB80qFyCmSPa0/m68vz+43WMmLKCY6ezvI4lQcDX22obA+2cc+cAzOxlYAHQnpxbbUUkiCQUjWb8gBaMm7eVZ2dtYMP3R3hlQAstYCiX5esIoxSQd52BeKB0boGcLvBUIuJ3ZsawTrV5a2ArDp04S089HS5X4Gth/A1YaWavmdlkYAXwrJnFA7P9FU5E/K9tUlk+frQ9dSsWZ/g7y/njx+s4e07P58p/8/kuKTOrBLTM/XKpc26331Llk+6SEsm/M1nZ/OmTdby+eActa5ZmdL9mlC8e63UsKQQFsh+GmV2X+5/NgUpAeu6firnHRCRExERF8LuejfjHvU1YlXGIbi8uJHX7Aa9jSQC57AjDzMY75wab2ZzcQ//xYufcTf4M5ysz6w50T0pKGpSWluZ1HJGgt/77Iwx9axm7Dp7k6W4NeKBNdd16G8IKase9V82sonOus3OuM/A6cAxYA/QpgJwFwjn3kXNucEJCgtdRREJC/Uo5t952qFuO38xcy+PTvtPe4XLFwhgHnAEwsw7AX8gpjcPAeP9GExEvJRSN5tUHkvnZLXX5cOUu7np5EekHTngdSzx0pcKIdM6dv4h5LzDeOfeec+5pIMm/0UTEaxERxmO31GHSgzeQcfAE3V5ayLxNmV7HEo9csTDM7PzDfTcDX+f5nq8P/YlIkOt8XXk+GtmeSgmxPPTat4yZs5lQWodOfHOlwpgCzDOzGcBJcp7uxsySyLksJSJhonqZeN7/SVu6N67Ms7M2MvStZRw9ddbrWFKILlsYzrk/AU8Ak4H27t9/pYgARvo3mogEmriYKEbd15Sn7qjP7PV76TUmhc17j3kdSwrJFZ/0ds4tcc594Jw7nufYJufccv9GE5FAZGY8cmMt3hzYkkMnztJrTApfrtNufuHA16VBRET+Q9vaZZk5sj01y8Yz6I1Unv9yE9nazS+kqTBEJN+qlCzKP4e2oU+LRF78Ko1Bb6RyRPMaIUuFISLXJDY6kmf7NOb3PRsyb1MmvUansHnvUa9jiR+oMETkmpkZD7SpwduPtOLIqZyl0met/cHrWFLAVBgiUmBa1SrDzBHtqV2+GEPeXMY/NK8RUgKmMMyslplNNLPpeY7VN7NxZjbdzIZ5mU9EfFO5ZFGmDWnDXc0TGfVVGoPf1PMaocKvhWFmk8xsr5mtueB4FzPbaGabzeyXAM65rc65gXlf55xb75wbCtwDtPNnVhEpOLHRkfz97sb8tnsD5mzMeV5ja6ae1wh2/h5hTAa65D1gZpHAGKAr0ADoa2YNLvUBZtYD+AT41H8xRaSgmRkPtavJWwNbcfDEWXqOSWHOhr1ex5Jr4NfCcM7NBy7cgaUlsDl3RHEGmAr0vMxnzHTOdQXu919SEfGXNrXLMHNEO6qWiuPHry/VOlRBzIs5jCrk7Np3XgZQxczKmNk4oJmZPQlgZp3M7EUze4VLjDDMbLCZpZpZamamVtEUCUSJpeJ4b1hbuuWuQzViygpOnMnyOpZcpYBZcdY5tx8YesGxucDcK7xvPLl7cyQnJ+uvLSIBqmhMJC/e15SGlUvwzOcb2Jp5nPEDWlC1dJzX0cRHXowwdgFV83ydmHtMREKcmTG0Y21eeyhnf40eoxeyeMt+r2OJj7wojKVAHTOraWYxwH3AzGv5QDPrbmbjDx/WiusiwaBTvfLMHNGe0vEx9J/4DW8s3q55jSDg79tqpwCLgXpmlmFmA51zWcAIYBawHpjmnFt7LT9He3qLBJ+aZeP5cHg7OtUtx69nrOXJ91dzOkv7hgcyC6VWT05OdqmpqV7HEJGrkJ3teO7LjYyZs4UW1Usxrn8LyhUv4nWssGJmy5xzyVd6XcA86X0tdElKJHhFRBj/c9t1jO7XjLW7D9Nj9EJWZ+jf5UAUEoWhS1Iiwa9b48q8N6wtEWb0GbeImd/t9jqSXCAkCkNEQkPDygnMGNGOJokleXTKCv72+QYtXhhAQqIwdElKJHSULVaEtx5pRd+W1Rg7dwuD3kjV4oUBIiQKQ5ekREJLTFQEf76zEX/o2ZC5mzK5c+witu877nWssBcShSEiocfMGNCmBm8ObMm+Y6fpOSaFlM37vI4V1lQYIhLQ2tYuy8zh7alQoggPTPqWySnb9JCfR0KiMDSHIRLaqpWJ4/2ftKNzvfL89qN1/OqD1ZzJyvY6VtgJicLQHIZI6CtWJIrxA1owvHNtpnybTv9Xv2H/sdNexworIVEYIhIezj/kN+q+pnyXcYgeo1NY//0Rr2OFDRWGiASdnk2rMG1IG7Kys7nr5UXMWvuD15HCQkgUhuYwRMJPk6olmTmiPXXKF2PIm8sY/XWaJsP9LCQKQ3MYIuGpQolY3h3Shp5NK/P3Lzbx2NSVnDqrFW/9JWB23BMRyY/Y6EheuLcpdSsU59lZG9mx/zjjH0imQolYr6OFnJAYYYhIeDMzhndOYvyAFqTtPUaP0QtZlXHI61ghR4UhIiHj1oYVeW9YW6IiIrh73GI+0oq3BUqFISIhpX6lEswY0Y7GiQmMnLKC577YqBVvC0hIFIbukhKRvM6veHt3i0Re+nozw99ZzokzWV7HCnohURi6S0pELlQkKpK/9WnM/91en8/X/sDd4xaz+9BJr2MFtZAoDBGRizEzBnWoxaQHb2DH/hP0HJPCip0HvY4VtFQYIhLyOl9Xnvd/0pbY6AjuHb+EGSt3eR0pKKkwRCQs1K1QnBnD29M0sSSPTV2pyfB8UGGISNgoHR/DW4+04p5kTYbnhwpDRMJKTFQEz9zVmKfu+Pdk+PeHNRnui5AoDN1WKyJXw8x45MZaTHwwOWcyfHQK36XryfArCYnC0G21IpIfN11XgfeGtSUmKoJ7XlnMx6v0ZPjlhERhiIjkV72KxZkxPOfJ8BHvrOCF2Zu0TPolqDBEJOyVyX0y/K7mibwwO41HtUz6RWl5cxERcp4M//vdjalToRjPfL6BnfuPM+GBZMprmfR/0QhDRCSXmTG0Y23G9W/Bpj3H6DkmhbW7dTPNeSoMEZEL3NawItOHtQHg7nGL+UJ7hgMqDBGRi2pYOYEZw9tRp0Jxhry1jJfnbgn7yXAVhojIJZQvEcu7g1tzx/WVeObzDfzP9FWcycr2OpZnNOktInIZsdGRvNS3GbXLFWPUV2nsPHCCcf1bUDo+xutohS4kRhh60ltE/MnM+NmP6jLqvqasTD/EnWNT2Lz3mNexCl1IFIae9BaRwtCzaRWmDm7N8dNZ3Dk2hYVp+7yOVKhCojBERApL82ql+HB4OyonFOXB177lzSU7vI5UaFQYIiJXKbFUHNOHtaFDnbI8/eEafvfRWs6Fwd4aKgwRkXwoHhvNqw/ewI/b1eS1lO088vpSjp4663Usv1JhiIjkU2SE8evuDfhjr0bMT9tHn5cXk3HwhNex/EaFISJyjfq3rs7rD7dk9+GT9BqTwvKdB72O5BcqDBGRAtC+Tlk++Elb4mKiuG/8EmZ+F3p7a6gwREQKSFL54nw4vB1NEhN4dMoKRs1OC6nlRFQYIiIFqHR8DG890orezarwj9mb+Nm7obO3hpYGEREpYEWiInnunibULl+MZ2dtJP3gScYPaEGZYkW8jnZNNMIQEfEDM2N45yTG9GvOml2H6TU2hbQ9R72OdU1UGCIifnRH40q8O6QNJ89k03vsIhakZXodKd9UGCIifta0aklmjGhHlVJFeei1pbwVpMuJBFRhmFktM5toZtPzHOtlZhPM7F0zu9XLfCIi+VWlZFGmD2tLx7rleOrDNfz+o3VBt5yI3wvDzCaZ2V4zW3PB8S5mttHMNpvZLwGcc1udcwPzvs4596FzbhAwFLjX33lFRPylWJEoJjyQzMPtajApZRuD30jl2Oksr2P5rDBGGJOBLnkPmFkkMAboCjQA+ppZgyt8zlO57xERCVqREcZvujfkDz0bMmfjXu4et5jdh056Hcsnfi8M59x84MAFh1sCm3NHFGeAqUDPi73fcjwDfOacW+7ftCIihWNAmxpMeugG0g+coNeYFFZnBP4GcF7NYVQB0vN8nQFUMbMyZjYOaGZmT+Z+byRwC9DHzIZe+EFmNtjMUs0sNTMzeO8+EJHw06leed4b1pboyAjufmURn6/5wetIlxVQk97Ouf3OuaHOudrOub/kHnvROdci9/i4i7xnvHMu2TmXXK5cucIPLSJyDepVzFlO5LqKJRj29jJembclYJcT8aowdgFV83ydmHssX7Snt4gEs3LFizB1cGtuv74Sf/lsA0++v5qz57K9jvVfvCqMpUAdM6tpZjHAfcDM/H6Y9vQWkWAXGx3JS/c1Y+RNSUxdms6Dk77l8InA2pCpMG6rnQIsBuqZWYaZDXTOZQEjgFnAemCac26tv7OIiASyiAjjiVvr8dzdTVi6/QC9X05h5/7A2ZDJAvVa2dUws+5A96SkpEFpaWlexxERuWZLtu5n6FvLiDBjwgMtaFG9tN9+lpktc84lX+l1ATXpnV+6JCUioaZ1rTJ88JN2JBSNpu+Eb5ixMt/TvAUmJApDRCQU1Swbz/vD2tK0akkem7rS8w2ZVBgiIgGsVHwMbw5sSe/mORsyPT7tO05nebPhta/IAAAFGklEQVQhU0hsoJRnDsPrKCIiBa5IVCTP3d2EWmXj+fsXm9h18CSvDGhBqfiYQs0REiMMzWGISKgzM0bcVIcX+zZjZcYh7hybwtbMY4WaISQKQ0QkXPRoUpkpg1pz9FQWd45dxJKt+wvtZ4dEYehJbxEJJy2ql+LD4e0oV7wIAyZ+w3vLMgrl54ZEYeiSlIiEm6ql43hvWFta1izNE//8jue+2Oj3O6hCojBERMJRQtFoJj/cknuTq7Iy/ZDfd/ALibukRETCVXRkBH+963pOZ2UTFenfMYBGGCIiQc7MiI2O9PvPCYnC0KS3iIj/hURhaNJbRMT/QqIwRETE/1QYIiLiExWGiIj4RIUhIiI+CYnC0F1SIiL+FxJbtJ5nZpnAjtwvE4CrbRBf3+PL6670mkt9/3LvKwvsu2I6b+Tnn3dhfn6wng+X+14gnw/g33MikM8HX1+bn//PL/e9azkfqjvnyl3xVc65kPwDjPfXe3x53ZVec6nvX+59QKrX/1wL8p93YX5+sJ4Pl/teIJ8P/j4nAvl8KIhzIlDPh5C4JHUJH/nxPb687kqvudT385M7EPg797V+frCeD1eTI9D4M3cgnw++vja//597dj6E1CWpUGdmqc65ZK9zSGDQ+SB5Fcb5EMojjFA03usAElB0Pkhefj8fNMIQERGfaIQhIiI+UWGIiIhPVBgiIuITFUYIMLNaZjbRzKZ7nUW8YWbxZva6mU0ws/u9ziPe88fvBRWGx8xskpntNbM1FxzvYmYbzWyzmf3ycp/hnNvqnBvo36RS2K7y3OgNTHfODQJ6FHpYKRRXc0744/eCCsN7k4EueQ+YWSQwBugKNAD6mlkDM7vezD6+4E/5wo8shWQyPp4bQCKQnvuyc4WYUQrXZHw/JwpclD8+VHznnJtvZjUuONwS2Oyc2wpgZlOBns65vwDdCjeheOVqzg0gg5zSWIn+IhiyrvKcWFfQP18nVmCqwr//tgg5vwyqXOrFZlbGzMYBzczsSX+HE09d6tx4H7jLzF4meJcSkfy56Dnhj98LGmGEAOfcfmCo1znEO86548DDXueQwOGP3wsaYQSmXUDVPF8n5h4T0bkhFyq0c0KFEZiWAnXMrKaZxQD3ATM9ziSBQeeGXKjQzgkVhsfMbAqwGKhnZhlmNtA5lwWMAGYB64Fpzrm1XuaUwqdzQy7k9TmhxQdFRMQnGmGIiIhPVBgiIuITFYaIiPhEhSEiIj5RYYiIiE9UGCIi4hMVhkghMrOfmlmc1zlE8kPPYYgUIjPbDiQ75/Z5nUXkammEIXKNzGyoma3M/bPNzOaY2a1mttjMlpvZP82smJk9ClQG5pjZnNz3vmxmqWa21sx+5+3/EpHL0whDpICYWTTwNTAeGAR0dc4dN7NfAEWcc7+/cIRhZqWdcwdyN8H5CnjUObfKo/8JIpel5c1FCs4ocgrjIDk7n6WYGUAMOev/XMw9ZjaYnH8XK+W+T4UhAUmFIVIAzOwhoDo5i8DdAXzpnOt7hffUBH4O3OCcO2hmk4FYP0cVyTfNYYhcIzNrQc4v/v7OuWxgCdDOzJJyvx9vZnVzX34UKJ7730sAx4HDZlaBnD2ZRQKWRhgi124EUJqcyWyAVOAhYIqZFcl9zVPAJnLmNz43s93Ouc5mtgLYQM4WmymFHVzkamjSW0REfKJLUiIi4hMVhoiI+ESFISIiPlFhiIiIT1QYIiLiExWGiIj4RIUhIiI+UWGIiIhP/h/j+m7EPKEzcwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# two methods can be used in the object\n",
    "\n",
    "# first is to calc equation 36 in Umetsu+15\n",
    "zeta       = 10.0**np.linspace(log10(0.05), log10(10.0), 100)\n",
    "sigma      = halo.Sigma(zeta = zeta) # unit should be Msun/Mpc**2. Since h = 1, so it is equal to Msun/h/(Mpc/h)**2\n",
    "# plot\n",
    "pyplt.plot(zeta, sigma)\n",
    "pyplt.loglog()\n",
    "pyplt.xlabel(\"zeta\")\n",
    "pyplt.ylabel(\"Sigma\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "triNFW/triNFW.py:216: RuntimeWarning: divide by zero encountered in arctanh\n",
      "  case_smaller_than_one   =   1.0 / (1.0 - X**2) * ( -1.0 + 2.0/np.sqrt(1.0-X**2) * np.arctanh( np.sqrt((1.0 - X)/(1.0 + X)) ) )\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'y [Mpc/h]')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAEKCAYAAAAM4tCNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJztnXvsbUd137/r9/PjGgwmPBTzUiAKokEur7gElD6S8KhJaChRo0IfFFB7lT8SiJqGNDEECEnFI0JNIFLjxk5JIaaogIzA4WECTSNiYwPGGAwJIEUQSEhIeRgb2/f+Vv84e/aZWbNm7Zk5+/z2OdfrI12dc/Ze8zjn3jvznbXWzCZmhuM4TisHS3fAcZz9xAcPx3G68MHDcZwufPBwHKcLHzwcx+nCBw/HcbpYbPAgohNE9BEi+gQRfYqIXrFUXxzHaYeWyvMgIgJwT2a+lYjOBvAnAF7EzNcu0iHHcZo4a6mGeTVq3Tp8PHv44xlrjrMnLDZ4AAARHQL4KIDvA/DbzHydYnMSwEkAOMThD9wD9z7eTo4dmbxgmJZty4UrypgmhZs9ZapMNu1vC0ZFs7WxRcwpcqb5c7Ia3eD2U9/Enadvr/oVF1u2JJ0gug+AdwD4WWa+uWR3b7ov/yA9eZsdUa4dDC+kXk8YbEjWc3BQrv8grWcsK9vT2gw2Rr9lPVnf4vKiL3q98rtZfZC/WfnfJGe/rzVAdN4rof3WtRx1/P+x/s813COrbVmP9Tmq58N/9Qf4xh1/XfWDLKo8Asz8dSL6IIBLABQHj62R/SM/iN4W7in/4HoGjbxM/UAw1tsykFkDWJVtqQ/Tg0c2QJTKVdan0jIQ9Aw0WXsNtuE/aUu78X/y7PdMTZPBJNgy25+B9XdoHAiXjLY8YFAcIKLzADwVwGeW6o/jOG0sqTweCOCNg9/jAMBbmfldx9qDwsydqY3YpmYZUJrJo2VBUXEYKqKokGpUUGlpEt+rUkGlvuS2VUuRhqVNUVXUzORzqAyr3prlvzVVt6iSgpqIf2+S9ZUUSHytUUosGW25CcDjlmrfcZzN8AxTx3G62AmH6bEz53LFWgaUIikxpeVK3F7PcqXFeVuK8Fi2Pc7QlujNVDnretyvbS1XBBu3UprGj/qqC87UyeVLfK0RVx6O43Rx91Ee8Ug7pTji+3Mqjp4wrNLfYhjW6IMZUi06ZMvqx3SGbhLy1ZBtV4Vsp03mhI+m+9SkToIaiH+nDifopAKJrzX20pWH4zhdnPnKw0ik6koAs7Iw51QcNeqnJ+GrM/xaVBpz+TFa1IU15R2TjyPjMHpf8CFY6qSq1/J7B3+I5ceYUiCxbSOuPBzH6eLMVR7HpTgU30emOOZO/Nok4asigtIUOen1oQxkCqNHVTSojd7oC82xB+xAqWNQD6FfVJPMJf0hVkRmSoHE9TbiysNxnC7OPOUxd8r5RKr5xrkbm0RSenI2aiIoDfkdLbtq1VlfTl81PhSrvqn6O2FUzPIDTfO4UCPSL2LWlSkQQ00I5aSmsjfiysNxnC588HAcp4szZ9myrZTzOZcrLYlfNWHYHmdlS/h1LmeoNkXta6g2hGQNB2opJFsXji0vYzJn6mhkJJSNhcvLmPDvgI5qOzk0V2/qOI6zZr+VhzHLbJRybsz65slfc5z41aI4ahK+rLTyqeSwhoSyRCmUnKA7EKo9FrSQLHRFUuz5qCby8pMKBMjDuKUyslwDrjwcx+liP5WHGs6bOKi4JgGsR3EYiiZTHHMlfm2S8KWFamvCryWlYfkzWkK1lnrI6imbruudQY1YIcwxqauhPkWRtKmJyjKyXGxj+WoafzNXHo7jdLF/yqMQsZhVcbScNWr6BToUR03il+UfkYqjwrYmglJUGjX+jK5oi1KmMDNu7cCfQyXSIdvW1MmUKjHUxFiF1a9BwVQpkNLJ6DP4Plx5OI7Txf4pD6Ds34jubT3lvEFxNKWa1+RWWLP9lG1FBKXKN1ETbanKBZHtVERkWoRGjSqpmHHHPkjTQzH7x2XEbE9Wx4M/RG6Us/pZ8oEkjQr1J6MvG+DKw3GcLnzwcByni/1bttDBfKecBzZZrmgO3JpHPU4tV6xkrpJTVGvLcopOhV97bUWZGudn1dKksASpCt1aiHrN8GswlbtUlQ6TcFaGfvYsX+Km142KPmi2NcuTToezKw/HcbrYL+VBwkk60ynns6WcTykObcNdyZFpJHNtFIbtDL9OKo3EVtR3YDhiw6WWRLCWmbJ07qfB2LYxa+fqQXGYDjZyV0LsSB1VTvZvMFUgaeMNtrKMdJwCnp7uOM7xspjyIKKHAvh9AN+N1bB9GTP/Zl3heU45nzvlvOlpbQ0JWpOKI/ludn1mSLUjVDvO0laI1fBfZArDartUfwuH5Vt5qDOqX8zqjJow7ODrEPfif15BhUwqEBgh2fHfkmI7XpjhDFbBksuWUwB+npk/RkT3AvBRIno/M396wT45jlPJYoMHM38FwFeG998iolsAPBiAMXhQvoAMd6YUx5ZSzklTE2O9RkRlym8R+xBaEr9KKqJiA1tNtKXKn1Hjx8h+q/Sj6R8xmFIj1inoxUQwYFQs6/KDbYiKKIXWaqQcmQn/nDIFMhrE/2ZEWvpYYYeqSA4OShqsr6K91fkhoocBeByA65btieM4tSwebSGi8wG8DcDPMfM3lfsnAZwEgBO4x+qilfa9LjjYWvkdMygOdY1ecXiP5bfAxKE9FbkbRRWxaQRFKo2CykjKGOqn5BcxFUSHq6OqXiEq1JTzzH8RIh3K9x8qHBXIuIEtry9TILJMDTUqZWx4z9PTiehsrAaONzPz2zUbZr6MmS9m5ovPphPH20HHcYosNnjQauq+HMAtzPy6pfrhOE4fSy5bfgjAvwXwSSK6cbj2y8x8tVnqgApp5GKpUHro9FDH6uIGKec9y5WaxC+1DxP11SxFmkK14b5Sb4szNEhxw/lZk5ZeTBjb9LSw7IyL4YVDu8n6QtiIMGwUAh6XO8OSobh8ieqbXL7EzOk43YAloy1/go1Wr47jLMniDtMWCFJBKOneJcVRSmtPjfP7U87ZmpRzTa1MJX5p/SrYNiV+xV99QnFY4VfLGVrlBM3upW2r90rUJItps/JhWm5UHOFCnEYuqwtKwUxLn1AgpX7tCTsRqnUcZ//YK+UxYp33iXBLztblZ6ZUpZyvL5T7MKU4tpX41R2qRdKW5dcoKQ1NZZT8GImCyLbkp0XM0G8PSn0yIUtuiIs3YUo1MppqW/FFWnpJgaxs1zWlZYarR5rtwFQ4ttpmePUnxjmOcxzsl/IgKid9yc1ulie+QnGobUftNKWcaxGV4078kr6KqK2S4khsDwu2st3oXq5sogJWtEapv3itlUhtZG2OOVxKf4UaYaFA0mpSf0hJgazuST+IdLzsLq48HMfpYr+UR0CJrBSjIZqq0KIqSfWaH6M/5Vw9LvC4czekcgA2iqDkn9e2mdLQlMyE0ujdGFdE5Gkk1WbH+YkyiNQIhAIJBkpkZlQcHRGV8fR0bft+xUFEx4ErD8dxuvDBw3GcLvZv2RI5TZtOO0/k+oFaRk09nyPl3EgSmz3xa+jvZBg2qWf4bDhFi8sUwxm6XraIV/n94nqt6WyOZUtcnTyhPDhM1d20oQviu3HZcbpuKNQhHKdAXQr7FMojK3uWSq248nAcp4v9Ux7FUG1BcZh1FRTHtlLOK8K6dWnkw+fDvN4pxaGGiwsJX7FtUWlYztCCc3VVN1KkrXKvC+OMDj5MbUYlgry/40Y4FmXC5/i3KiSSydBtfC1DqBXTVi2/faeqKw/HcbrYP+WBSCnMtdmtFI5N7jVMf9tKOZdKw0r8OijYKL4ULqmUWCkcluoTr1o9B9M26+vImdPXkXwpUb/YGUdHsWmqHtZ9S30WAMqJZNrkP6aGV/g6ZL0tYiK0c7qhTGWVjuM4TeyX8qBBLWyYem76OCb7UPZjFA/02TTlfEpxHCh9mIqkAJmfYVQIVgTlQC+r+jMsm8wvgpRS0l8rnM7Saco50nvSNxG7vo4KNuNXi/4ORCIZTgdbRb0UNudVRV2owkYw9mHfzzB1HGd/2S/lIdGe7GYojun6FIWwyWY34/pkyvlhud4wO/Fh8NHkbVXlbkiFIfwlWgRlrWzC5/S+Wl5TF+I3ynwoCqUt+ebaf/yxwuf8CMBStCVWCOF7jn4Q4S+Ju1Da6q/1U1UjMYpSym1owmA7uPJwHKcLHzwcx+liz5Ytg7PU2ilbKqmGX493p2zL7lfduYrBRixXNKetSPySy424rSzhywitjolVYhmjLnFkvaYNElpODVOXDDIMO16nzHYMsx4I2zjdWzpR2VhuiPDtWMu4xIn7INqeMZS6bVx5OI7TxZ4pD6CYDCbv79Bmt6aUc9W5OrwWFEfitJ1K/DpU6s2coJpKCWV0h6kV1tXURVlxNEgODRmClEqE83uUbU7L+5Y5UYWKSNLepYNUhGxNjPDrpHP1mHHl4ThOF/unPAB7m/0ObnarOuW8lP6NCsURnzYlFEeW+KVudkNaRvg10nrDZ9mXyFaoCl15kG6D3HaSRE1IKRNeB1VxlN/CGIZNFUiyMU4UGkO3p6VBVLFQEWrIVvpk9ghXHo7jdLFfyoMwRFvCEK4kiW2t7ZKqMPwjph9DVyWZfwOYVBypmrAVR7qJLi3Psr5EeQjbTIEotpa6MFRJZivIkq0UP8ZoK7fZRwfxrO8NCAUSNzS+k34HTTmIa5z/nJNYiWWbMOchQa48HMfpYlHlQURXAHgGgK8y80Uz1RnerF7nerLbVNnofVlVxP6R8KofG6imnBcUByf9TeuRqkL1Y4hIylplaH1I72URm6iNpmhLhfLIniWrpXvLQ3qk7VH8d5D6QdbqIu8TyTdZCvraeHJ2V/wjuxZJqWFp5fE/AFyycB8cx+lg0cGDmf8YwN8t2QfHcfrYeYcpEZ0EcBIAThycv5LPLWd07MBOWStJrLSj1Uw5L4VjlXqyszVi5+phek06QRNbuZQR9Wmh2tw26ufUMkVZtsjVQJaKHl3jEJIVu2GTPgjHaLY0if/OswdDpQ5YUtLKt+X07OJo2qSVpZctkzDzZcx8MTNffM7BiaW74zjOwM4rj5RBdchEsO7qCoqjZrPbQM1mNy1JbK0mhDrRNpHJJK5CODaut6gQDlG0zT/H39O2SdWPLJP2Jf2+SLH+SrNUc+lBXb8nYTJWH8/A0jEqTxujsq2pJgpO1a0xV/iVuanPO688HMfZTZYO1V4J4IcB3J+IvgTgZcx8eWdlq9eWbfZTdWnlS6pCuVc6ezSuj0X98gSw2CZL3tJSzicUh+7HmFYTUzaprWhLmcnlc1tqtuBnPo5QKFYTYhZePxUu70Pm2wi+j56JPK53BiEwm5+kRpUc9TW26ODBzM9Zsn3HcfrZM58HgAPa/jZ70Z5axtpEV4riGGv+/OAcxT8SbA/TfmtRjJLiqFETR4eabfqKDpWizfo1Z5cWt9UHxRGXHZPAhAKRdcUXp16BTJ2Uktt6qUkbz2ysMlP1+enpjuMsxf4pj156t9nLMpqvY2Bqs1vqdxH9yvwkyGzlsYFq2veE4tB8Hkcl2yQyA/WermiGN4VU9qTP8ruEy0oEJUs9D2W0HAaZZFETQVGqOePwjXGO4yyNDx6O43SxX8sWGhLDZFgWqHOUFustlNXKlxLBlHvF1PPYRp4GJhyecT3ZDlZtKTKxXImXF1PLlWQpIq6FJVSoA4ZzVU0kCxT+ejR1PSZ4iWVKUsV4bqj4bC1bgulx7Ww1H1JVOIM1ohTGtZyuc57jEXDl4ThOF/ulPIDVjD3XGR2ZqugJvyp9yBK/lD6MqgSJDUvVgjw0u3Yy1oRf089HWnp6pirSz6tywrZQJrlHqU2iTiYEYZJGLidNWTa+nykOoSYqhGjTyV/GqexVzCEIrE1vNRviOPZI13eoOHgQ0eMryt/FzJ+sbs1xnDMGS3n8HwDXwx6AHw7gYXN2aHZqttmXttujHJqt2WZvn+EplEYh9TztJxJbNeVc2mR+hzxUK9XD6AuJU+QPSzZpHVp9a99M/mwTybiej0/8kslgYfu7VBlx+UIoWG22J0Zb8wiW0qlmGlIhaD4KmUbekyQWrnempMdYg8f1zPyjVmEi+qONe+A4zl5SHDymBo5am3kZ/B0tB/yY1VVEZlq22UvVQMJGjcykr5l/I653tEFio876wt+gppEX1USqbFSbs0S9qi2n94xIx1h2mIEpmhnptAyLDNelykCjv6KAqRQ2mLDXSsSKisjPxkGttden7nVS5TAlogcD+J7YfjhC0HGcuymTgwcRvRrAvwTwaaxXmwxgtwaPUjTE2hhXKKs+L3Ui6hJTPNhHi4pkfpL0NS5X2uI+12a3zJcSqYmjs9JrUoEcnRX5M0qRmANrZhy+0mlNOoSt8sNFsc1e++vqipiYNrZRohDkYUW15Sb61JTfcWTcy9oS+f+V1CiPfw7gkcx8R1PNjuOc0dQkiX0BwNnb7ojjOPuFlefxeqzE020AbiSiDwAY1Qczv3D73VM7pqeR15ST70vnkarh11LINvqQJXEV6ojfZ07WXG/nD2MSS5AkAU6+irKqY1NfrhxF/zrkskUuVzixTR2lPDhO1akqW64MF05RZhNCvRufXVsgPy9Ds0ltqSYJKzh/W9LSZVn1XmUdal+2G6q9YXj9KIB3btyS4zhnFNbg8UMA/hDANcz8rWPqTx2bng5WKJNdt+4p6kU9XQzIw7LR+6ozOkrPV1EcpqXnqZhnjU4kgCU2QXGcvZq5jsRnAOCzgxoJikNRHjIZbLyXOkUBZOeQbhSO1U5al5Nwg7NSS0+fTA6LrxfUgxmi7TlRrMZ525g4Zvk8LgfwGABXE9EHiOgXiegxTbU7jnPGYiWJXQfgOgAvJ6L7AXgagJ8nor8P4OMA3sPMbz2ebg4Q6pK/gFwhaAllA5ucDpZujENyL0tBVxLKarbZl87WlEpkdU1/lad6AZHqkQokqIlYeZQUxzmDyjgnVh7DdDq8HpyVOwZ4UBx8atU4hw6GEGP8BLYOH4c8ab0qRVzYaqeZqRvhKvtChqooh2GtesXNI+NeQFMinf6PqiQxZv4agCuHPyCiH4A/oNpx7tbUJIn9FwCvYeavD5+/C8CzmPkl2+6cygHpakLbACfLCdtJf4Z2T6uvVL6Qgp6UzzbPKclcUu3IqE7L81XizW4Fn8d4XYmgSMVxdO4w3Z2znvYOzl3JhrPOGV7PWr1yJJFO3bVq7FToy2njtx8o+hK0jXEyFdzwTZR9H5pCSKMsWf1xuZ4oy+gLMnweLVGWGTfCSWryPJ4eBo5VX/j/Afix2XviOM5eUbNsOSSic0OGKRGdB+Dc7XZrAyrSxwOT2+3jexPb7eN7VQf8lKIrmkoRCsN8mn12T38FpqMsY34Gcp/H6OMYFMfheWsnxYnz7gQAnH9ilRZ01sHK5o5T639u38Y5AIDTo88jNCRek2uUfM6eIGfcs/wYJXUS2xbVjmFbLKOe9l7whxjqJ3s6Xo0fQylLce5Hg0CpGTzeDOADRPR7w+fnA3hjfROO45yJTA4ezPxqIvoEgKdgNS69kpnfu/WeOY6z09SeYXojVkuVIwCfmKtxIroEwG8COATwu8z8qokSaXq6di6ppGNXrXpP1qfdz+oVr9bpYDIFXUk5R2Yj6k9s08/m2R+lE8Vih+kYxh2WLSEMOzhHw1IFAO5//rcBAA+65zdWZYfOfOlb9xltwrKFxVIkW6Igdk6mnyE/a/eMEOvk0sYK1TaQhWjjJUOhXjVEW0qFN88wrU8Sa2XSYUpE/x7ARwA8C8C/AHAtEb1g04aJ6BDAbwN4OoBHAXgOET1q03odxzkeapTHLwB43JDrgSFh7MMArtiw7ScA+Bwzf2Go9y0AnonVuSH9zOEo7T2XVNRTOqsjeV/hMC05V2uSxLJwrqE8xjytcAJYbDsqjmEGG5RHCMcG5yiwVhxvedJlAIBn/+lJAMCpo3WFp08PjtIQoh2SxsIGuSRJTCaOCcWhOkFlmFRRKWV1wul1VIRo45m90LblKM3Klu7H1yzn6oTDNElF7zzPoyZU+zUA8d6Wbw3XNuXBAL4Yff7ScC2BiE4S0Q1EdMOdp2+boVnHceagRnl8DsB1RHQVVmPoMwHcRET/EQCY+XVb7B+Y+TIAlwHABedeyLHPQz0RXVLjz7A2xgUsX8doE+6JV9Xvsnop+Tq0k9bz9HckZbW282e9rE0zdSJPO082xg2zUtgWP3w+PBwUyMF6Wg0+jqA4vvztCwAAt35nHeE/dedQ+V2rxuiuQXEMW/FjhXAwqhMk97QnyJX8IZpKkeqkyedhhICn0tGtcPG67TwhrJSObp4W1qkqaqgZPD4//AlcNbzea8O2/xLAQ6PPDxmuOY6zB9SEal+xpbavB/AIIno4VoPGswH8q41rrTkoqLRYUzfGyfrFdeWAH/UM1Cl6fB5WQllBnejb90U7oyLhzDa8kggHxAlgIaoSfBxBcXzn9nNGm6M7VsqDgvIYFMfB4Nc4ODWajr6OTEVIBYK1mihGZDSFMN4r+DW0a5lCiH0IeVuJbe8zWaZ8HZrfRUBa9GULT4wzDwBi5p+obkUvf4qIfgbAe7EK1V7BzJ/apE7HcY4PS3k8CSuH5pVYbc3vmE5tmPlqAFfPXa/jONvHGjwuBPBUAM/BajnxbgBXLqoOCABR3eMUxPWqxykYdU2eEhb6l5TRXxPbiTKarRUunlzamOFiLtpm/RxCq2F3bEj6AoDbaHVe9qlTw87ZwTkalioAgDuH5cqdw3JldJgOzUUPegrLlrCkIfmqOExLjtLkYVITyWeWw1Q6StWwrnB6ao7YKUdp92MV5FLmOJPEmPk0M7+Hmf8dgCdiFXX50LDUcBznbo7pMCWicwH8OFbq42EAfgvAO7bfrQ1oObt0LFO4ntQ3UUdyr9CXZNbXQ7RW/W0OUzu8q9Vbevj06ubwKk4AC37NsDs25ihcG8Oxa5tRcYTXoaKDU6kCASKFUFAcdFpRE9mrcKRG32m8J8OmWpi05CjVVIpwxHY5Sq3EL1mvcWZH5ijV6m0888NymP4+gIuw8km8gplvbqrZcZwzGkt5/BsA3wbwIgAvjJ6VQQCYme+95b6VMU4Sk1SdT5pdLzct72nPYjFncEnBJ2GmshsU1UkLQmUA69l+TCBDml6etDMmb4nEr7vWRsHHMSoO4fM4iNLT136Q4VUojkRNTIRzrSSxPJU9Dn2m13KfR6R+xOxePN0subeBr2My/FrJXOd5MHNN6rrjOHdTarfk7wzFrfeBio1xLVvym57FkpVNX63ohSyT9q9g06BILOQsyuPMGxmNSmBoLNzTfudRGQxqYkz8ipRMwcch/RrAdJTFjrak6iL1Y4jXoCrEBjzV9nRab9OzWGo2pVnKoUZVyA1wVvSl83zTorogoo9NFa6xcRznzMRSHt9PRDcZ9wnABTP3pw3jWSym7VimYGvV1RJlMTA3wsV1xWW29IzW9aw8qIAwE8XNBdUQZrTTRl+OUsUhoyVArjikulDT08d8D1avJ22J6IpMQTdtG05EVze/FXwdEP6ShKzt1AeStZH0pcO/McOWfGvw+HsV5U9PmziOcyZiOUz/4jg74jjOfrF3DtPkDFOD2UK0RcdmxRKicE5pVRnTZtpkNJUKXHECjmpVODoT41Gml36Q/P16mZEuY+L341KhsHxJbfWlSHr2R8mmnCSWhWa1sG4hNCvLJOVqHqdgnG+alW1wrlY5SjfEw7GO43RR87jJnwXwpuFJcbvP1IliFfc0hVCccJWHV3fRkJjWVb028cjwpYxgRw3zcIJYqS9p8lUazjU3sE1selu9Z9PGDNWyfn31Xg/NqpvoSqFZxREL6XCVjtIkTCrL5Eomq7flRPRSCruann4UdXqaGuXx3QCuJ6K3EtElRNty+TuOs0/UnCT2EiJ6KYCnYfW0uDcQ0VsBXM7Mn7dLb5mWM0xr7lkPrx5tprtVtdFsrL/CJmvAqE5OLNm6Xmk7vMrkqNh2+G1K3U3Cj6XZ3/J5FMKxqk0pEcy4p6sJcU0qA+V0sHV96Q+9tXNJa5LErBPRt0iVz4OZGcBfDX9OAfguAP+biF6zxb45jrPD1Pg8XgTguQD+FsDvAvgFZr6LiA4A/DmAF2+3i0lvqqMt5rA4h5t4k81vjeWziIlEubGeuWQ6eWQk19XG70LaGrzUdsVZo6HttR9DREm0xK9MeaS+EO2eFW2RyWBrH0iuJqYP+OHMtpTEZUVQMsVhpZGXUtC18lZC2dFRfq2CmlDtfQH8pMz7YOYjInpGU2uO45wx1Pg8Xmbcu2Xe7mxOOYejPNXXpIaXfSp1/ZqDMCvxmDiC9FW7JhVHPPPKBozoi/yNap4BK5+rokdF0oiHGhWZUBzW6el5/Xmko7jNXjsRPdtwp/gbStEVKzJT2EzXnUY+pTiOjjJbbtyS73kejuN04YOH4zhd7F96eom7UfqJtWTIT8UKMlXbppu+kgwxx5HwQlnr0YwtIdXioySj91Op59q1YjgWyjLFOKNjfU++hv7HDlOkNi2ng9WcYTp1Vscx4crDcZwu9lt5bLrRrCItXTLbmRol9SBDrKateAXWM1mpwbjeg/RWy8lk8kyKmg1nVc9XqUr8suq1FYduK23KIdXSeaRpklh9+LV4Tzvdq+d0sJKjNLJlq00DVx6O43Sxn8qj4qSuuXwgc2xKW1e2WfFsLc0iZBu9DyZZclc8RY6KIz06wPzp5HcQM3ByTSZSVSmPfCafUhw1KeeabXYqmFQccYq88G0UE8Gi90U/SaLSCopjk7CsWm9fIpjFIsqDiH6KiD5FREdEdPESfXAcZzOWUh43A/hJAL/TVIowj6I47siM9Cko1ERQRhVk2EqlwamogJoinzkwlA5Kk4ICSfpQUiBQlECmEGJboTSkUjid2xYVh/qsWqkmlO+Y2RT8D1H5kuLQoy7Tfow5oyuc9DfU17Ylf5HBI2Sm+u5+x9lfdt7nQUQnAZwEgBNnVTykbodcwOOGNm2n/1FqlEVxkjV06tsgaRPlcIRT2cf6wyntYdarOHVA7W9JabByX1x2/ZqRAAASqElEQVTLnj6P3LdRd8p5qDdVHFoaeUlxaBvjsuMHFTWR1SP9T+qGu4Li0HI35Gdrm/2UnyS5l/o61MgK9/lDtjZ4ENE1AC5Ubl3KzFfV1sPMlwG4DAAuOO+B83l7HMfZiK0NHsz8lG3V7TjO8uz8smXXGHe2khImnXruVJyYk60V0lcy7o3+13FpEju/pIc0/ciaw1TcUzcVF5YtWnp6MZFKSSQrOkFVW+E4NU4uLy1XSF2KiM8y9TyuJ1wLywHNCTq1XLGSuWTfZjqjgxuXJDUsFap9FhF9CcCTALybiN67RD8cx+lnqWjLOwC8Y9Y6txS5sZye1WVGqRAZZUojLZTME+PMnT6LJlMgwKhCwpkfeZJY4X3+MUVMXCV1kV4rfI77JVSJdZbG6CA1QqpTisM6z0MqDjq9/vEmN7upqewTDk7t3rY3tyVhWdGHI26J1O5SbMJxnH3CfR69KGoifzpbmMLK/pEwM/JhWq/m82DYCiSubz2Fi35vqNDyp5+lzWnX1PBrdi+tP078yp7oJn0gsfKQm9wK4di0Pl1xqJvdxvNNLX8D9HsVfgwzEWzK11E6HSwpo4RlGzfEBVx5OI7ThSuPiCySEjOO1BUzdxjJDwe/g+I3KZ1HGnwfHLcjRESmQNSUc9GnUip6JX3RFp60ydSE9GvE96SfZCyj2QqFoG2Mm1IcTX6M+LtNKA6t3m1vs2fp/FKosYlw5eE4TheuPFoJg/Pgo0hyN5AqgbUPRJEeCKZhhhiiLAfienRvVA+jryNf88vt9ZFcST4WupPXJ9F8HIDqz8h9H5aaEDaWSjGeApdFf+QzWDRFU1Ic8Xea2F6v53lMq4lJxdGwzb4ql8NKTz9itIRbXHk4jtOFDx6O43Rx91y2yBBqdj96X5T2Zedq6fGQ2iFe2bIiKNFoWJdO1OyhTMlaRKxPjNBsU9BWpk8XlibpNeGQ1FLOSzbaMqiUap4scezlStXDmcZwsbVkkH1QbAtlN94pW0MpNNviQJ3AlYfjOF2cMcojC33WMKVAIhsa6s2cjPGsFxRCqC/MbAch5TyaTUN9LMoqaeXSiTrahOqUx0Ku7zXMVhKlaDlkW+EMVWzykG8+65c32qUqIClfUhzW81Vk26oTVHyXmpBqaK9ms1sNPWd0WIlgR1E9Dd1w5eE4ThdnjPIYEaHUGkYfhUjuUhlHciq3I9wOLBTIyiQ0mralJYmNqdpjqFY2EFfQstttAlV56A6dGn9GVSKZenq68IvIxCw19Clt83p7Us63vtmtZ5u9tdlNkoRlg2mfOnXl4ThOF2ee8gjI0bRhQ5iV+FVjm/kmwv14dhKp61KIUDL1ivrlCelxZEZGW3omFaOMdh5p0m58r+DPSK5lEY9yfU2JXyXFEW8em1QTyvdriYpMpZ5b9chEsOha0xPepCoxfR+enu44zjFw5ioPQaIQLH9FCRl1EWnlAEDiIJ6657QMZcYiymFAY5lUVYSoTrJ9fTQutz2J5eso2TT4M1IbUb+mPAob4uTGtrhczdPaphSHmnI+pRTifh7XM2WtbfYFxZE+t6Xv9HRXHo7jdOGDh+M4Xdxtli0qMoknLBlieR12sHacwCUT15IaTovlikgkoyTzC4mtrD+9mJbpQiubOUjl/Xpn6OqavgwyHzJdWtpUJH7VLUWETYtjUw0XTzhZ1XsNO2WtRLBC+jnLdkTfW3Dl4ThOF3dv5dECp28yxykwOk+l41Q9HWygmEgWO8qEM1WqgESRbKI4Qp9U5cG2jeEMXV+PZ0ZxreQUjdsuKA7dtkVNFGxqHJvSOWrYtDwWUmVqs5tRRlUcme1Rkwpx5eE4Thf7pTwYq5FRPMQ5Zj2ThxPG85F0yib2N6wVQehC6lRIbEX4tqRA0q8k/CHCF7KqJ1Ujue8jq7YuKW5illHr1WZ7oM2fAXSFX6e20Mf1zp74NX6XguLoDesWFEfXZrfYzyEVR+l7xH1oxJWH4zhd7Jfy2EViL3shgSxLHoso+UPMyIys5EBRGS0e9MKaWfd9SJtcKWT3rAN+IO6pKmVCcRxj4tfYzlwp54KtbXaTCkfpAzM3ucxceTiO08UiyoOIXgvgnwG4E8DnATyfmb9+bB0I69Xw0fSPpLkfme9DOQKwqEDicV3UY21oG9PQwzUpNE63zBcGRjVVfoyAFTlBeq8pgiI/L5G7UWqnJjKjKQ6Zz2E92S2rP1Ucqn+jUnEk9VaylPJ4P4CLmPnRAP4MwC8t1A/HcTpZZPBg5vcx86nh47UAHrJEPxzH6WcXHKYvAPC/Nq5FOrTCZW1J0pJqLpcixdBt3Ghh+RIjnKljc5pT1VjSZN2d+G7qMiOrpKKcVL8NztCkPiGVrUc9ViVzif7NnviVfSfle0wtV6K65liuqEwknbH2OzSytcGDiK4BcKFy61JmvmqwuRTAKQBvNuo5CeAkAJw4695b6KnjOD1sbfBg5qdY94noeQCeAeDJbOz+YebLAFwGABecuJC1hCAAerhyupOrvoSPFSpFnndK0cIvVw2pAknvDBWN55MKp6q2oAxtW89iaQnRTpXRIopS4WlKYbQ16m8Iv2b9bAnVyn5XpJGbiV+l+lSlpCsOddYvKY7iYyHjapX2CqHa6Y129f9+loq2XALgxQD+CTPftkQfHMfZjKV8Hm8AcC6A99NqFr2WmX961hZKqgKITh+XqoLLtqNRGjbVtu+PXTCcFFKNZH6R+HQw2U9rdpCKpSXz2JiVavwYa9sJf4bWVkv4NdCT+KV+h4LisJTHtlLOs/rzlPO8PstPUlAcFdv3p1hk8GDm71uiXcdx5mMXoi39aP6PMPNaUYfRVlcVLbbJ82flTKtFZNbWqa116GhQRDVuHe080wLFh8kZ/ot1WWMmk/Vbs76mHmSZLltDGQjbprNGZTsV/oamlPOKiMosikNRG3zELS4PT093HKeP/VMe8UjasO08ySKX5Qo5IoltSYEk5Ya8jlCBEpEZ65VKQ+ZyaPXLpmfYdi86JdrTZjTDj1Fq04qKSJueaIvWv5rcjVLbBZ+NblOe/efO4Zg80Ce63qw4CvcsXHk4jtOFDx6O43Sxf8uWElLuGo7TLCQrh9BI8o3B1o7ly1id5kwck8Hk8kW0Z9GREJb0oZgcZjg/S22rUtxopyf8Kj9veuJXqS81tjWnnM+wXFEfzlToS1UfQtEZQrWuPBzH6WI/lYc2wsqZOoys2vA4hlmlArHUhOiCploKM/n6HA4lTb10ildy9scmz44cqNj8NKkutHuWY1PatqgTK/xq1Sud3w0h26rw7hyKozflPOvvBopDUxuNStaVh+M4Xeyd8qAjBlsz8RjytMOxsY3pA5HlSmViDAWTdcFQJSNznRQmKc001gwkJqwef4ZariX8apRpOvGrVH/FrG+evjVnyrnSdk/K+ZyKI+DKw3GcLvZOeSRYvo8s+hI/+7VBTRwUyoSiWrc0NQI0KZKt0aAqYkyfRKnemmiL/NySxBX6Zm2z76m/R3FoKeclm5YEMCXa0pUAljeoX2/AlYfjOF3sl/JgAMygsKSsUBMj8QisPA82LhPPlJmKqFAPJe9FUZGs7g71T1bfRkUIv5jvUaMYWqIt2rWaCIq4V/RrWOU7cjdSkwk1kSiEGRXHpinnVj83xJWH4zhd+ODhOE4X+7VsEcTyNQvfWsuYMYGssHyJi4nPPcuYUl0qw3kcVenpWhstsnTboVr5ucZZafShK+Frqh2gL/Frqr3YtuWhTOJ6d8q57MOMy5WAKw/HcbrYP+URj6DR7BxmJZYb4mpmUalASHHEynCu3CgXM0NItucU9CqONVTb4DitCb9OlEnaLKqfijCsrAtoS/wS94oPnVbr233FEXDl4ThOF/unPGIUFVKlQEqJZIGKIdUMuhZDsnG4eLqNWekJ2VqzljVLl2y0kGpouyb8Kj+3+FA2CcOq9W2gOAzfR1fKeU1/t4ArD8dxutgz5cFFn0fmmygpEMU2+1yTUGZQEyexE8bmo+7ZtA1RlxY/hmEzuYHNKt+S8CWud0VSYpvxXkXi11h2C6eca305BrUR48rDcZwu9kx5CCw/RkmBtFLyFWw47G5dd2yiOKyDg1pSzwfME8vl5xa105A3sa3cDbX+ilTzYr/2QHEEXHk4jtOFDx6O43Sxf8uWI9ZPRp9wgqqp7KUyFppTVVIzJGsJab20yNaKs0yrnJbGdTPsKq9tEvK1TtsaP1YsQTrCsHn90yd/ZWUq2jRPOV9ouRJYRHkQ0SuJ6CYiupGI3kdED1qiH47j9LOU8ngtM78UAIjohQB+BcBPV5fWTkYvhV2VGV7OjJkzVVM0NdSoEskmw3dLO4GWxC/jnuoELdXTs4HNulaT8DXaGiqiJjlM1GOmmmf1FhSHpYLGjxX9XphFlAczfzP6eE+g5dncjuPsAqSuv46jYaJfB/BcAN8A8CPM/DcFu5MATg4fLwJw8/H0sIr7A/jbpTsRsWv9AXavT94fm0cy871qDLc2eBDRNQAuVG5dysxXRXa/BOAEM7+sos4bmPniGbu5Ed6faXatT94fm5b+bM3nwcxPqTR9M4CrAUwOHo7j7A5LRVseEX18JoDPLNEPx3H6WSra8ioieiRWyd9/gfpIy2Xb61IX3p9pdq1P3h+b6v4s5jB1HGe/8fR0x3G68MHDcZwu9m7w2LXUdiJ6LRF9ZujTO4joPgv356eI6FNEdEREi4UAiegSIvosEX2OiP7zUv2I+nMFEX2ViHYiT4iIHkpEHySiTw9/Xy9auD8niOgjRPSJoT+vmCzEzHv1B8C9o/cvBPDfFu7P0wCcNbx/NYBXL9yf7wfwSAAfAnDxQn04BPB5AN8L4BwAnwDwqIV/l38M4PEAbl6yH1F/Hgjg8cP7ewH4syV/I6yOmDl/eH82gOsAPNEqs3fKg3cstZ2Z38fMp4aP1wJ4yML9uYWZP7tkHwA8AcDnmPkLzHwngLdgFZJfDGb+YwB/t2QfYpj5K8z8seH9twDcAuDBC/aHmfnW4ePZwx/z/9beDR7AKrWdiL4I4F9jtaluV3gBgD9cuhM7wIMBfDH6/CUs+B9j1yGihwF4HFaz/ZL9OCSiGwF8FcD7mdnsz04OHkR0DRHdrPx5JgAw86XM/FCsslN/Zun+DDaXAjg19Gnx/jj7ARGdD+BtAH5OqOpjh5lPM/NjsVLPTyCiiyz7nTwMiHcstX2qP0T0PADPAPBkHhaNS/ZnB/hLAA+NPj9kuOZEENHZWA0cb2bmty/dnwAzf52IPgjgEhgbUXdSeVjsWmo7EV0C4MUAfoKZb1uyLzvE9QAeQUQPJ6JzADwbwDsX7tNOQUQE4HIAtzDz63agPw8IkUIiOg/AUzHxf2vvMkyJ6G1YRRPG1HZmXmxWI6LPATgXwNeGS9cyc/3BRvP351kAXg/gAQC+DuBGZv6nC/TjxwD8V6wiL1cw868fdx9Ef64E8MNYbYH/awAvY+bLF+zPPwTwfwF8Eusz+n+Zma9eqD+PBvBGrP6+DgC8lZl/1Syzb4OH4zi7wd4tWxzH2Q188HAcpwsfPBzH6cIHD8dxuvDBw3GcLnzwcCYhoocR0e1D6nK4xkT0pujzWUT0N0T0rg3aeSIR/Xcieh4RvaFg80EiunXJHcPOCh88nFo+P6QuB74N4KIhoQhYJRVtmm/zdADvsQyY+UcA3LBhO84M+OBxN4eI/sFwFskJIrrncJaDuach4moAPz68fw6AK6N6X05E/5OI/pSI/pyI/kN07xeJ6JPD2RGviup7MoBrhvcPIqL3DGVfs8FXdLbETu5tcY4PZr6eiN4J4NcAnAfgTcxce2DOWwD8yrBUeTSAKwD8o+j+owE8EaujEz5ORO8G8BisthX8IDPfRkT3BQAiuj+Au5j5G6vMbTwWq52mdwD4LBG9npnjnbrOwvjg4QDAr2K1H+U7WB2wVAUz3zRsJ38OVipEchUz3w7g9mGj1ROwGlx+L+wDYuZwxsbTALwvKvsBZv4GABDRpwF8D9Jt/s7C+LLFAYD7ATgfqxOtTjSWfSeA30C0ZImQex+svRDS33FH9P40fKLbOXzwcADgdwC8FKsjDl7dWPYKAK9g5k8q9545+FLuh9WmtOsBvB/A84noHgBARPcddpg+GsCNSh3OjuKj+d0cInouVr6GPyCiQwAfJqIfZeY/qinPzF8C8FuF2zcB+CBWO1lfycxfBvBlInosgBuI6E6sljtvB/Dx4zgLxZkP31XrTDL4Nd7FzLVRGBDRywHcysy/UWH7EqzOPH1LZd0fAvCfmNlDtgviysOp4TSAC4joRpHrMQvM/Gu1toPj9XsB3DV3P5w2XHk4jtOFO0wdx+nCBw/HcbrwwcNxnC588HAcpwsfPBzH6eL/A/T6bbmfsuvdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Second is to calculate the projected map\n",
    "\n",
    "# make coordinates\n",
    "xmpc_edges        = np.linspace(-3.0, 3.0, 100) # Mpc/h\n",
    "ympc_edges        = np.linspace(-3.0, 3.0, 100) # Mpc/h\n",
    "xmpc_bins         = 0.5 * (xmpc_edges[1:] + xmpc_edges[:-1])\n",
    "ympc_bins         = 0.5 * (ympc_edges[1:] + ympc_edges[:-1])\n",
    "xmpc_mesh, ympc_mesh  =  np.meshgrid(xmpc_bins, ympc_bins)\n",
    "\n",
    "# mapping\n",
    "map_is            = halo.Sigma_XY(XX = xmpc_mesh, YY = ympc_mesh) # ordering is (ny, nx)\n",
    "\n",
    "# plotting\n",
    "pyplt.imshow(map_is[::-1], \n",
    "             interpolation = \"nearest\", \n",
    "             norm = LogNorm(vmin=1E12, vmax=1E16), \n",
    "             extent = [xmpc_edges[0], xmpc_edges[-1], ympc_edges[0], ympc_edges[-1]]\n",
    "            )\n",
    "pyplt.xlabel(\"x [Mpc/h]\")\n",
    "pyplt.ylabel(\"y [Mpc/h]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}