small_scale_dust_pysm3_generate_templates.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "papermill": {
     "duration": 0.239907,
     "end_time": "2023-02-13T04:22:37.550932",
     "exception": false,
     "start_time": "2023-02-13T04:22:37.311025",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "os.environ[\n",
    "    \"OMP_NUM_THREADS\"\n",
    "] = \"200\"  # for jupyter.nersc.gov otherwise the notebook only uses 2 cores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import healpy as hp\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_dir = Path(\"production-data\") / \"dust_gnilc\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "datadir = output_dir / \"raw\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "output_nside = 8192"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_lmax = int(min(2.5 * output_nside, 8192 * 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Large scales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/global/homes/z/zonca/condanamaster2/lib/python3.8/site-packages/healpy/fitsfunc.py:643: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  alm.real[i] = almr\n",
      "/global/homes/z/zonca/condanamaster2/lib/python3.8/site-packages/healpy/fitsfunc.py:644: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  alm.imag[i] = almi\n"
     ]
    }
   ],
   "source": [
    "alm_log_pol_tens_large_scale = hp.read_alm(\n",
    "    datadir\n",
    "    / \"gnilc_dust_largescale_template_logpoltens_alm_nside2048_lmax1024_complex64.fits.gz\",\n",
    "    hdu=(1, 2, 3),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_log_pol_tens_large_scale = hp.alm2map(\n",
    "    alm_log_pol_tens_large_scale.astype(np.complex128), nside=output_nside\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.73506452,  3.735317  ,  3.73569709, ...,  3.75179248,\n",
       "         3.75241809,  3.75209597],\n",
       "       [ 0.01206824, -0.01207586,  0.01203825, ...,  0.02188258,\n",
       "        -0.0218166 ,  0.02191317],\n",
       "       [-0.01173623,  0.01173273, -0.01174043, ..., -0.02359359,\n",
       "         0.02363389, -0.02364048]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "map_log_pol_tens_large_scale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Small scales modulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "modulate_alm = {\n",
    "    k: hp.read_alm(datadir / f\"gnilc_dust_{k}_modulation_alms_lmax768.fits.gz\").astype(\n",
    "        np.complex128\n",
    "    )\n",
    "    for k in [\"temperature\", \"polarization\"]\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Small scales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl_small_scale = hp.read_cl(\n",
    "    datadir / \"gnilc_dust_small_scales_logpoltens_cl_lmax16384.fits.gz\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/global/homes/z/zonca/condanamaster2/lib/python3.8/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return array(a, dtype, copy=False, order=order)\n"
     ]
    }
   ],
   "source": [
    "synalm_lmax = 8192 * 2  # it needs to be the same for all output nside\n",
    "# synalm_lmax = 1024\n",
    "np.random.seed(8192)\n",
    "\n",
    "alm_log_pol_tens_small_scale = hp.synalm(\n",
    "    list(cl_small_scale),\n",
    "    lmax=synalm_lmax,\n",
    "    new=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "alm_log_pol_tens_small_scale = [\n",
    "    hp.almxfl(each, np.ones(output_lmax))\n",
    "    for each in alm_log_pol_tens_small_scale\n",
    "]\n",
    "map_log_pol_tens_small_scale = hp.alm2map(\n",
    "    alm_log_pol_tens_small_scale, nside=output_nside\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16385"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmJklEQVR4nO3deXjU5bn/8fedkBAIYQsJAQIMIAliURFE1Fp3BWPFWmvBrVZOPXpqV7sEq9bjRk6X03O0tlaLpZ5fq7XUWmuiuIsLLVuVRQiGEEhYwxYSwpLl+f2RBDMxCRNmMt9ZPq/rmuti7tlunivJZ77b85hzDhERiU8JXjcgIiLeUQiIiMQxhYCISBxTCIiIxDGFgIhIHFMIiIjEsR5eN9AVgwYNcj6fz+s2RESiyvLly3c55zLaeyyqQsDn87Fs2TKv2xARiSpmtqmjx7Q7SEQkjikERETimEJARCSOKQREROKYQkBEJI4pBERE4phCQKSLdu4/xI79h7xuQyQkFAIiXeCcY8pDr/PZ/3oD5xzle2rRmhwSzaLqYjERrz3yRgkAdQ2OUXOKADg3J4OHrppAUoKR2TfFy/ZEukwhINIF//3q+k/V3l5fydkFbwBw34yTuO6MkSQmWLhbEzkuCgGRAO1sPg5w2YQsfnXdJKoP1fHSqu384C8rjz7nnr+tYeu+Q3zn4rH07JHoVasiAVMIiAToo237AbjprFEApKUkcc3pwxncL4WSnTV8UL6Pv3+4lcfe3sBjb2/g6knZ/OxLp3jZssgxKQREAlSyswYA36DefvVzczI4N6dpgsabzvLxxV+/D8CC5RUsWF7B0H4pFH7zHAakJoe3YZEA6OwgkQC9uHIbI9N7k9GnZ4fPmTRyAH/7+tmMzexztLa16hAT73+V3TWHw9GmSJcoBEQC4Jzj4x3VnJ+biVnnB31PGd6fV797LvO/ejr9eycdrU964DUu/cUi1u+o7u52RQLmaQiYWYKZPWhmj5jZV7zsRaQzZbtrOXCkAV9672M/udl5uZl8cM8lPHPL1KO14h3VXPKLRUx+4DWqautoaNQ1BuKt4w4BM3vSzHaa2eo29WlmVmxmJWaWf4y3mQEMA+qAiuPtRaS7vb52BwAXjR/c5ddOHZ1OyYPTeXjWxKO1XTWHOeW+VxhzZxGH6hpC1qdIVwWzJTAfmNa6YGaJwKPAdGA8MMvMxpvZBDN7sc0tE8gFFjvnvgvcFkQvIt2qYu9B0nr2YFj/Xsf1+h6JCVxxylA2zr2M71yU4/fYuLtf5p2PKxUG4onjPjvIObfIzHxtylOAEudcKYCZPQPMcM7NBS5v+x5mVgEcab6r3wCJWBV7DzJsQK9jHg84FjPjWxeNZeaU4fzhH5t4uPkK5BvmLQEgZ3AfXvzGOST30OE6CY9Q/6QNA8pb3a9ornXkOeBSM3sEWNTeE8zsFjNbZmbLKisrQ9epSBds2XfwuLcC2jO4bwrfvSSX+2ac5Fdfv6OGnLte4oPyfSH7LJHOhPo6gfa+JnV45Ms5VwvM7uwNnXOPA48DTJ48WUfRxBMVe2s53Tcg5O9745k+bjzTx9pt+5n+v+8crV/56HsALLvrIgZ1ckqqSLBCvSVQAQxvdT8b2BrizxAJq70HjlB9qJ7hAwI/M6irThzSlw/uuZgXbj/brz75gdfw5ReytvlqZZFQC3UILAXGmtkoM0sGZgIvhPgzRMJq5ZYqAE4a1rdbP6d/72ROzu5PWUEev71xst9j0//3HXz5hdQ1NHZrDxJ/gjlF9GlgMZBrZhVmNts5Vw/cDiwE1gLPOufWhKZVEW+UVjZNF5E7OC1sn3nR+MGUFeQxZ/o4v/rYH73E9/78oa4vkJA57hBwzs1yzg1xziU557Kdc/Oa60XOuRzn3Bjn3IOha1XEGxV7D5KSlMBAD+b++fdzx7Bx7mVc0ur6hAXLKxhzZ5GuPJaQ0HloIsdQvqeWYf2DPz30eJkZj9842e/KY4BLfrEIX36hwkCCohAQOYa12/czLqt7jwcEYurodD6452JuOsvnV28Jg5UV+zzpS6KbQkCkE4fqGijfc5DcrPAdD+hM/97J3HvFSZQ8OJ3pn8nye+yKX76HL79QVx5LlygERDpRWd00/XNWhK0d3CMxgV9fP4l1909j9mdH+T027u6X+cvyCg4eURjIsSkERDqxszkEMvpG5gVbKUmJ3H35eP5198X0Tfnk2s87/vwhJ97zMr99p9TD7iQaKAREOrGjeV3hzhaSiQQDUpNZee+lfPzgdL/6A4Vr8eUX8u7HuzzqTCKdQkCkE8vK9tKzRwIntFopLJIlJSZQVpDHH//tDL/69fP+iS+/kNLKGpzTNQbyCYWASCeKd+znxCF9SUlK9LqVLjnrhEGUFeTx5vfO86tf8PO3GTWniPdLtGUgTRQCIp3Y0jyFdLQaNSiVsoI8vn3RWL/6tb9t2jKoPlTnUWcSKRQCIh1obHRs3XeI7BBOIe2Vb1+UQ1lBHk+0mZNowr2vMOn+V6k9Uu9RZ+I1hYBIB3bVHOZIQ2NUbwm0dfH4wZS0OXi8+8ARxt+zkCUb9+h4QRxSCIh0oGLfQQCG9oudEICmawzKCvJYee8lfvVrfrOYUXOKKFq1zaPOxAsKAZEOlO06AED2wNgKgRZ9U5IoK8jj6+eP8av/xx9W4MsvZMXmvR51JuGkEBDpwNvrKxmYmszYzMiYMqK7fP/ScXz84PRPLXV51a/ex5dfyK6awx51JuGgEBBpR0Oj4811O7lgXCaJCd7MHhpOSYkJ3Himj3X3TyNvwhC/x1pWN2vUGgYxSSEg0o5dNYfZf6ieU7L7ed1KWKUkJfLodae1u5to9J1F+PILtaBNjFEIiLSjZeK4jLTImjgunL5/6ThKH7rsU/UxdxbxxCLNSRQrFAIi7dhZ3TRnUGaEThwXLgkJRllBHuvun+ZXf7CoaU6it4p3etSZhIqnIWBmI8zsBTN70szyvexFpLWWLYHMtPgOgRYpSYmUFeSx4u6L/eo3/W4pvvxCNu+u9agzCVYwC80/aWY7zWx1m/o0Mys2s5IA/rDnAIXOuZuB8cfbi0io7dzfFAKDInz20HAbmJpMWUEef7ntLL/65376Jr78QmoO68rjaBPMlsB8wG8b0cwSgUeB6TT9UZ9lZuPNbIKZvdjmlgn8C5hpZm8AbwbRi0hIVdYcpl+vpKibOC5cJo0cQFlBHtdPHeFX/8yPF3Lm3Nd15XEUOe4QcM4tAva0KU8BSpxzpc65I8AzwAzn3Crn3OVtbjuBrwI/ds5dAOQdby8ioVTX0MhTizdR19DodSsR74ErJ7Bx7mWcm5NxtLat6hCj5hRx/W//qTCIAqE+JjAMKG91v6K51pGXgW+a2WNAWXtPMLNbzGyZmS2rrKwMWaMiHSmtbLpSeEi/+D0zqCvMjN/fPIUND11G60sq3i3Zxag5Rfx04TrvmpNjCnUItHdVTYdfBZxzq51zVzvnbnXOfa+D5zzunJvsnJuckZHR3lNEQqp8T9NBzp996RSPO4kuiQlG6dw8Pvyx/5xEj765AV9+IRsqazzqTDoT6hCoAIa3up8NbA3xZ4h0q83NITB8YG+PO4lO/Xo1zUn0ync+51e/8Odv48svZHvVIY86k/aEOgSWAmPNbJSZJQMzgRdC/Bki3Wr5pr0M7ZdCemqy161EtZzBaZQV5PHad/3DYOrc1/HlF7Kv9ohHnUlrwZwi+jSwGMg1swozm+2cqwduBxYCa4FnnXNrQtOqSHis3lrFxJEDMIv9OYPC4YTMpjC49/P+Z4Gfet+rXPWr9zzqSloEc3bQLOfcEOdcknMu2zk3r7le5JzLcc6Ncc49GLpWRbqfc47tVYcYqoPCIXfT2aPYONd/GooVm/fhyy9kznOrPOpKNG2ESCv7D9ZzuL6RwX0VAt3BrGkairZh8PSSzfjyC/ntO5qTKNwUAiKtbN/fMmeQQqA7dRQGDxQ2zUlUsrPao87ij0JApJWSnU2nMY5KT/W4k/jQEgbFD/hPUHfRfy/StNVhohAQaaV4+34SDMYO7uN1K3GlZ4+mCep+3ubajDHNaxhI91EIiLSybns1vkGpmjPII1+clE1ZQR4PfWGCX92XX4gvv1DTUHQDhYBIK8U7qhmXFdtrCkeDa88YQVlBHicO6etXHzWniBvm/dOjrmKTQkCkWe2RejbvqSV3cN9jP1nC4qVvnfOpg8fvfLwLX34hSza2nb9SjodCQKTZ+h01OAe52hKIKC0Hj1fd6z8n0TW/WYwvv5Cqg3UedRYbFAIizYq37wfQ7qAIlZbSNCfRwm/7T0Nxyn++gi+/kEN1DR51Ft0UAiLN1m6rpldSIiM0cVxEy81Ka3d1s3F3v6wziY6DQkCkWfH2anKy0khI0JxB0aBldbPbzhvjV285k0gCoxAQoWnOoOId1YwbrF1B0eaH08Z96uAxNIXBX5ZXeNBRdFEIiNC0pvCeA0d0UDhKtRw8Lnlwul/9jj9/iC+/kLJdBzzqLPIpBERo2hUEOigc7XokJrQbBuf97C2dSdQBhYAIn4SAtgRiQ0sYvJ9/gV+95Uyi+oZGjzqLPAoBEZqmixjUpyfpfXp63YqE0ND+vSgryON/vnyqX/2EH72kg8fNFAIiNG0JnDhEWwGx6sqJwygryOPM0el+dV9+ITMfX+xRV5FBISBxr6HRsX5HNbk6MyjmPX3LVMoK8vxq/yjdE9cL2igEJO6V7T7A4fpGHQ+II50taFO+p9ajrrwRthAws9FmNs/MFrSqpZrZ783sCTO7Lly9iLS2dlvLdBGaOC6etJxWuuLui/3q5/zkTXz5hRypj4+DxwGFgJk9aWY7zWx1m/o0Mys2sxIzy+/sPZxzpc652W3KVwELnHNfA67oUuciIbKsbC+9khIZp2MCcWlgajJlBXm884Pz/eo5d70UF6ubBbolMB/wW//NzBKBR4HpwHhglpmNN7MJZvZim1tmB++bDZQ3/1uzP4knlmzcw2kj+5OUqL2j8Wz4wN6UFeTx1M1T/OqxvrpZQD/1zrlFQNvJu6cAJc3f8I8AzwAznHOrnHOXt7nt7OCtK2gKgg57MbNbzGyZmS2rrKwMpF2RgFUdrGPt9v1M8aUf+8kSFz6Xk0FZQR4zTx/uV/flF3Lf3z/yqKvuE8xXn2F88i0emv6gD+voyWaWbmaPARPNbE5z+Tngi2b2a+Dv7b3OOfe4c26yc25yRkZGEO2KfFrJzqY1BCZk63iA+Cv44smUPuR/8PjJ9zbiyy9kaVnsLGjTI4jXtjfVYoc7z5xzu4Fb29QOAF8NogeRoOzYfwiAIf16edyJRKKEhKaDxw2NjjF3Fh2tf+mxpmsLlt91UdRfYBjMlkAF0Hp7KRvYGlw7IuG1raolBFI87kQiWWJzGHxwj/+ZRJMeeA1ffiGNUXzwOJgQWAqMNbNRZpYMzAReCE1bIuGxveogKUkJ9OuV5HUrEgX69246k+it753nVx8dxQePAz1F9GlgMZBrZhVmNts5Vw/cDiwE1gLPOufWdF+rIqFXvKMGX3oqZlpIRgLnG5RKWUEec6+a4F/PL+SCn73lTVPHKdCzg2Y554Y455Kcc9nOuXnN9SLnXI5zboxz7sHubVUktBobHf/avJfTRg7wuhWJUrOmjKCsII+JI/ofrZXuOoAvv5BfvVXiXWNdoBOjJW6VVNZQfaie00YoBCQ4f/2Psz81J9FPXi6OimkoFAISt1pO8zut1bc4kWC0t6BNyzQUtUfqPeqqcwoBiVvPrdjCqEGpjBqU6nUrEkNaFrRZe5/fJAuMv2dhRJ5JpBCQuNTQ6PiwfB+XnpSlg8LSLXolJ1JWkMcr3/mcXz3SziRSCEhc2lZ1kPpGx4iBvb1uRWJczuA0ygryyJ8+zq/uyy/kykff86irTygEJC6V7zkIwPCBulJYwuPWc8dQVpDn98Xjg/J9+PILeXZZeSev7F4KAYlLH1bsA8CXruMBEl6LfnD+p84k+sGClfjyC9lVczjs/SgEJC79dcUWTvcNYLh2B4lH2lvdbHLzNBThXMNAISBxp66hkdJdNUz2DfS6FYlzLaubrbvf/0yicK5hoBCQuFO+p5a6BseYjD5etyICQEpS05lEr99xrl/dl1/Y7WGgEJC4s6HyAABjMnQ8QCLLmIw+lBXk8b1LcvzqvvxCnl6yuVs+UyEgcae0sgaA0doSkAh1+wVjKSvIIy3lkyVf5jy3qls+SyEgcWdDZQ2D+vTU9NES8Vbde+nRM4ky07pn8ZpgVhYTiToNjY5lZXsZm6mtAIkeJw3t220LH2lLQOLKPzfupnTXAWZOGX7sJ4vEAYWAxJWPtu4H4LMnDPK4E5HIoBCQuLJ2WzWZaT2jfnFwkVBRCEjccM6xtGwPJ2f387oVkYgR1hAws9FmNs/MFrSqXWlmT5jZ38zsknD2I/Fl464DbN5Ty7m5mV63IhIxAg4BM3vSzHaa2eo29WlmVmxmJWaW39l7OOdKnXOz29Sed859DbgJ+HIXehfpkjeLKwE4LyfD405EIkdXThGdD/wSeKqlYGaJwKPAxUAFsNTMXgASgbltXn+zc25nJ+9/V/N7iXSLdz6uZExGqiaNE2kl4BBwzi0yM1+b8hSgxDlXCmBmzwAznHNzgcsDeV9rWtapAHjJObeincdvAW4BGDFiRKDtivhxrmklsUvGZ3ndikhECfaYwDCg9WoIFc21dplZupk9Bkw0sznN5W8AFwFXm9mtbV/jnHvcOTfZOTc5I0Ob8XJ8tuw7yN7aOj6jg8IifoK9Yri9xVk7nAjbObcbuLVN7WHg4SD7EOnUqooqAE4ephAQaS3YLYEKoPWll9nA1iDfUyTkVm2pokeCkZuV5nUrIhEl2BBYCow1s1FmlgzMBF4Ivi2R0Fq1pYrcrDRSkhK9bkUkonTlFNGngcVArplVmNls51w9cDuwEFgLPOucW9M9rYocH+ccq7ZUMUG7gkQ+pStnB83qoF4EFIWsI5EQq9h7kH21dUzQQWGRT9G0ERLzFm/YDcCpw/t724hIBFIISMwrWr2N4QN7MX5IX69bEYk4CgGJeWu27mfqqHSarksUkdYUAhLTDh5poLL6MCPTNVWESHsUAhLTyvfWAmi+IJEOKAQkppXtOgDAyPRUjzsRiUwKAYlpf1+5jbSePcgZrIXlRdqjEJCYdbi+gYVrtvOF04bROznYabJEYpNCQGLW6i1VHKlv5KwxWlRepCMKAYlJzjkee7sUgMm+AR53IxK5FAISk94qruTVj3bwzQvHMqhPT6/bEYlYCgGJSc8s3czgvj35xgUneN2KSERTCEjMaWh0LN6wm/NzM0lK1I+4SGf0GyIx56Ot+9l/qJ4zx6R73YpIxFMISMx5f8MuAIWASAAUAhJz3t+wm7GZfchMS/G6FZGIpxCQmHKkvpGlZXs4S1sBIgFRCEhMWVmxj9ojDdoVJBKgsIWAmY02s3lmtqBNPdXMlpvZ5eHqRWLX+xt2YwZnjFIIiAQioBAwsyfNbKeZrW5Tn2ZmxWZWYmb5nb2Hc67UOTe7nYd+CDwbeMsiHXt/wy7GD+nLgNRkr1sRiQqBbgnMB6a1LphZIvAoMB0YD8wys/FmNsHMXmxzy2zvTc3sIuAjYMdx/w9Emh2qa2DFpn06HiDSBQFNreicW2RmvjblKUCJc64UwMyeAWY45+YCge7aOR9IpSlEDppZkXOuMcDXivhZvaWKIw2NnO4b6HUrIlEjmGMCw4DyVvcrmmvtMrN0M3sMmGhmcwCccz9yzn0b+CPwRHsBYGa3mNkyM1tWWVkZRLsS6z4o3wfAqcP7e9qHSDQJZpL19lbtdh092Tm3G7i1g8fmd/K6x4HHASZPntzh+4ss2biHof1SyOyr6wNEAhXMlkAFMLzV/Wxga3DtiByfolXbeOWjHXz+lKFetyISVYIJgaXAWDMbZWbJwEzghdC0JRI45xyPvFHCCZl9+P6luV63IxJVAj1F9GlgMZBrZhVmNts5Vw/cDiwE1gLPOufWdF+rIu3bUFnD2m37+cpZPnpo1lCRLgn07KBZHdSLgKKQdiTSRf/cuAeAz56gZSRFukpfmyTqLSvby6A+PfGl9/a6FZGooxCQqFbf0Mii9ZWcOSYds/ZOWBORzigEJKotLt3N7gNHyJuQ5XUrIlFJISBR7cUPt5GanMh5ue3OTCIix6AQkKh1pL6Rl9ds5+Lxg0lJSvS6HZGopBCQqLVs0x6qDtZx2YQhXrciErUUAhK1lm7c27R2wGjNGipyvBQCEpXmvbuRX7y2nuwBvejXK8nrdkSilkJAotKv39oAwLSTdFaQSDAUAhJ1nHMcrmvgsglZ3HGJ5goSCYZCQKLO7gNHqD5cz+m+gTorSCRICgGJOht3HQDANyjV405Eop9CQKJOSwiMSlcIiARLISBRpbSyhiff3Uh6ajLZA3p53Y5I1AtmeUmRsNq5/xAX/PxtAH78+fFaO0AkBPRbJFGjZd2Am88exQ1TR3rcjUhsUAhI1FiycQ+pyYncedk4bQWIhIh+kyQqVB+qY9HHlZw2coACQCSEwvbbZGajzWyemS1oVUswswfN7BEz+0q4epHo850/fUjF3oNcr91AIiEV6ELzT5rZTjNb3aY+zcyKzazEzPI7ew/nXKlzbnab8gxgGFAHVHSlcYkfh+sbeOfjSm48cySXapoIkZAKdEtgPjCtdcHMEoFHgenAeGCWmY03swlm9mKbW0crfuQCi51z3wVuO77/gsS6VRVVHK5vZKpmCxUJuYBOEXXOLTIzX5vyFKDEOVcKYGbPADOcc3OBywP8/ArgSPO/GwJ8jcSZlrOCTvcN9LgTkdgTzDGBYUB5q/sVzbV2mVm6mT0GTDSzOc3l54BLzewRYFEHr7vFzJaZ2bLKysog2pVotWTjHnIG92FgarLXrYjEnGAuFrN2aq6jJzvndgO3tqnVAm2PE7R93ePA4wCTJ0/u8P0lNtU3NLJ8016unDjU61ZEYlIwIVABDG91PxvYGlw7Ip8orazhf177mJrD9UwZpeMBIt0hmBBYCow1s1HAFmAmcG1IuhIBfvTX1Swu3Q3AFB0PEOkWgZ4i+jSwGMg1swozm+2cqwduBxYCa4FnnXNruq9ViSe7ag6ztGwPyT0SmDVlOFn9UrxuSSQmBXp20KwO6kVAUUg7EgH+9sFW6hsdr3zrHHIGp3ndjkjM0vX3EpEWLK/g5Ox+CgCRbqYQkIizZmsVa7ft5+pJ2V63IhLztJ6ARIzfvbeRPy0t55yxg+iRYHz+ZJ0WKtLdFAISMd5Yt5N126vZuOsAZ45JZ4AuDhPpdtodJBHBOcfKiioADtc3cokmihMJC4WARITyPQepOlhH35SmjdOLTxzscUci8UG7gyQirNrStBXwyLWnkWDougCRMFEISERYuWUfSYnG1NED6dkj0et2ROKGdgeJ5xobHQtXb2fi8AEKAJEw05aAeGrR+krWbd9P2e5a7rgk1+t2ROKOQkA8s2ZrFTc+uQSAof1SmP4ZnREkEm4KAfHMT14upl+vJL53aS7jstLokai9kyLhphAQTyzftIe311dy52XjuGHqSK/bEYlb+uolnvjde2WkpfTgegWAiKcUAhJ226oO8tLq7cw8fTi9k7UxKuIlhYCE3W/eLqXROW6Y6vO6FZG4p69hEhZVtXXMe7eUIf178dTiMq4/YyQj0nt73ZZI3FMISLerOVzPV363hA/K9wGQkdaT70/TNQEikUAhIN2qsdFx+x9XsGpLFQ/PmohzjjEZfeibkuR1ayJCGEPAzEYDPwL6Oeeubq6NAH4J7ALWO+cKwtWPhMdvFpXyVnEl91/5Ga44RYvEiESagA4Mm9mTZrbTzFa3qU8zs2IzKzGz/M7ewzlX6pyb3aacAxQ6524Gxnepc4l4Rau28bNXisk7eQjXnzHC63ZEpB2Bnh00H5jWumBmicCjwHSa/oDPMrPxZjbBzF5sc8vs4H3/Bcw0szeAN4/vvyCR6Jklm/n6H1cwcXh/Cq6agJl53ZKItCOg3UHOuUVm5mtTngKUOOdKAczsGWCGc24ucHmAn/9V4MfN778A+F3bJ5jZLcAtACNG6NtkNFiztYq7nl/NOWMz+M31k+iVrJlBRSJVMNcJDAPKW92vaK61y8zSzewxYKKZzWkuvwx8s7le1t7rnHOPO+cmO+cmZ2RkBNGudLfNu2vZsu8gdzz7IQNSk3l45qkKAJEIF8yB4fa2711HT3bO7QZubVNbDVwdRA8SIX791gb+6+V1R+/P+8pk+vfWQvEikS6YEKgAhre6nw1sDa4diRZH6hv507JyThvRn6raOn66cB0Xjx/MebkZDE5L4UKtESwSFYIJgaXAWDMbBWwBZgLXhqSrKLG75jAvrtzGxeMHM7R/L6/bCasFyyu4+/mmk8V6JBijM/rwP18+ldSeuvREJJoEeoro08BiINfMKsxstnOuHrgdWAisBZ51zq3pvlYjT+Gqbfz4hTWc85M3ufX/lrN4w26c63CPWMxwzvHU4jLGZaVx34yTuGBcJr+5YZICQCQKBXp20KwO6kVAUUg7iiINjU1/8K87YwQvfLiVl9dsJ2dwH24808cXJg6L2T+KKzbvZd32auZeNYFZU0Zw45k+r1sSkeOkWUSD0PKl/46Lc/nHnAv5ydUnk5SYwF3Pr2bqQ6/zn39fQ2lljbdNBuHjHdX8cMFK3t+wy6/+1OJNpPXswYxTdQWwSLSLza+qHkhJSuSaycP50qRsVmzex1OLy/h//9jE794r49ycDL5y1kjOzckkMSHyL5pyzvH0knLue3ENh+qaDgBfOC6Ta04fjnNNVwJfd8ZIrQUgEgP0WxwKrf6umxmTRg5g0sgB/CjvRJ5ZUs4f/rmJm+cvY8TA3twwdSQzTh1KZt8U7/rtRPWhOn74l5UUrdrOOWMH8dAXJvDiym386s0SXl+3E4AEQyuCicQIhUAQjnUIODMthW9eOJbbzhvDwjXbeer9TTxYtJYHi9aSnppMblYauVlpjMtKI2dw0607jyNs3l3LP0p3s3b7foq3V7O96hCXTRjCTWf7GNSnJ+t3VHPr/y1n055a8qeP45ZzRpOQYNx23hiuPWMEm3YfoNFB35QejM7o0219ikj4KATCICkxgctPHsrlJw9l7bb9LN6wm+Lt1azbUc0zS8o5WNdw9LkjBvYmNyuNEzL7MKx/L7IHNN2G9e8d1NW3xdurmfHouxyqa6RXUiI5WWkM7d+LR98q4Yl3Spn2mSxe/WgHvZN78Md/O4MzRqf7vb5fryROzu5/3J8vIpFJIRCEltNBuzI32olD+nLikL5H7zc2Osr31rJuezXFzbd12/fzVvFO6hr8tzXSU5PJHtCLSSMHctt5Y8hI6xnQZ9Yeqefrf1xBWkoSf/v6GYzN7ENC87GJDZU1PLGolOdWbGFCdj9+dd1pDI7QXVUiEnoKAY8lJBgj01MZmZ7KpSdlHa03NDoqqw9TsbdpPp6KvS23Wn6/uIw/Ld3Mv587hn87Z9QxD9De/fwaNlTW8IfZZ5Cbleb32JiMPhR88WTuzDuR1OQeUXHgWkRCRyEQoRITjKx+KWT1S2Fym8c2VNbw05eL+e9X1/N//9jEdy7K4ZrJ2fRI/PQZvwuWV/CXFRV868KxnHXCoA4/Tyt9icQnhUAIhPu785iMPjx2wySWb9rD3KJ13PnXVcx7t5QvTspmwrB+TBjWj/69k/l4RzV3P7+aqaMH8s0Lx4a5SxGJBgqBKDZp5ED+fOuZvPLRDn7x6np+8nLx0ceyB/SirqGR3smJ/O/MidrNIyLtUggEoeWKYS9XzTIzLj0pi0tPymJf7RHWbN3Pqi1VrNpSxYadNdx9+Xgd6BWRDikEYkj/3smcfcIgzu5k37+ISGuaOygI7piXi4mIRDaFQAhob7uIRCuFgIhIHFMIBOGTA8Pe9iEicrwUAiIicUwhEAQdFhaRaKcQCAHToWERiVIKARGROGYt0yFHAzOrBDa1KvUDqto8rXWt7eODgF2EVns9hOI1nT2no8eONR5t74djfDrqK9jnh2N82t6PpvE51vMCGYv2atEwPoG+Jp7GZ6RzLqPdR5xzUXsDHu+s1vZxYFk4egjFazp7TkePHWs8vBif4xmjSBmfdsYrasbneMYoVsanu36GYml8Wt+ifXfQ349Ra+/xcPQQitd09pyOHjvWeLS9H47xOZ7PiZTxCbSXYHXH+BzreYGMRXu1aBifQF8Tz+NzVFTtDgqWmS1zzrWdnl+aaXw6p/HpnManc5E6PtG+JdBVj3vdQITT+HRO49M5jU/nInJ84mpLQERE/MXbloCIiLSiEBARiWMKARGROBa3IWBmqWb2ezN7wsyu87qfSGNmo81snpkt8LqXSGVmVzb//PzNzC7xup9IY2YnmtljZrbAzG7zup9I1Px3aLmZXe5VDzEVAmb2pJntNLPVberTzKzYzErMLL+5fBWwwDn3NeCKsDfrga6Mj3Ou1Dk325tOvdPFMXq++efnJuDLHrQbdl0cn7XOuVuBa4CIOzWyO3TxbxDAD4Fnw9ulv5gKAWA+MK11wcwSgUeB6cB4YJaZjQeygfLmpzWEsUcvzSfw8YlX8+n6GN3V/Hg8mE8XxsfMrgDeBV4Pb5uemU+A42NmFwEfATvC3WRrMRUCzrlFwJ425SlASfM32yPAM8AMoIKmIIAYG4eOdHF84lJXxsia/BfwknNuRbh79UJXf4accy84584C4mKXaxfH53xgKnAt8DUz8+TvUA8vPjTMhvHJN35o+uN/BvAw8EszyyN80ydEonbHx8zSgQeBiWY2xzk315PuIkNHP0PfAC4C+pnZCc65x7xoLgJ09DN0Hk27XXsCReFvK2K0Oz7OudsBzOwmYJdzrtGD3uIiBNqb7N855w4AXw13MxGoo/HZDdwa7mYiVEdj9DBNXybiXUfj8xbwVnhbiUjtjs/Rfzg3P3ytfFo87AapAIa3up8NbPWol0ik8Tk2jVHnND6di+jxiYcQWAqMNbNRZpYMzARe8LinSKLxOTaNUec0Pp2L6PGJqRAws6eBxUCumVWY2WznXD1wO7AQWAs865xb42WfXtH4HJvGqHMan85F4/hoAjkRkTgWU1sCIiLSNQoBEZE4phAQEYljCgERkTimEBARiWMKARGROKYQEBGJYwoBEZE4phAQEYlj/x/7VaK7ndKiMgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.loglog(hp.alm2cl(alm_log_pol_tens_small_scale[0].astype(complex)))\n",
    "len(hp.alm2cl(alm_log_pol_tens_small_scale[0].astype(complex)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_log_pol_tens_small_scale[0] *= hp.alm2map(\n",
    "    modulate_alm[\"temperature\"], output_nside\n",
    ")\n",
    "map_log_pol_tens_small_scale[1:] *= hp.alm2map(\n",
    "    modulate_alm[\"polarization\"], output_nside\n",
    ")\n",
    "assert np.isnan(map_log_pol_tens_small_scale).sum() == 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.11883597, -0.13409836, -0.13284604, ..., -0.03759012,\n",
       "        -0.04488804, -0.03506894],\n",
       "       [ 0.00683536, -0.00612315,  0.00449982, ..., -0.00329152,\n",
       "         0.00467237, -0.00257872],\n",
       "       [ 0.00955996, -0.00842392,  0.01196066, ...,  0.00195204,\n",
       "        -0.00244644,  0.00073119]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "map_log_pol_tens_small_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fac6975ebb0>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAn3UlEQVR4nO3de3hU1b3/8fc3gRDuCAQQCIRrEPCCpmC9ay0CoWqtbb3UY5VK1aJtbe0J6jm19UJ6fm3P0WqLWC3aU1HKsRYbVFRUUKkSQOUOIUQSAiSQcA0ht/X7IwlNxkmYZCazZzKf1/PkeTJrz575sJ4h39lr772WOecQEZHYFOd1ABER8Y6KgIhIDFMREBGJYSoCIiIxTEVARCSGqQiIiMSwDl4HaIm+ffu6lJQUr2OIiESV1atX73POJfnbFlVFICUlhezsbK9jiIhEFTP7vKltGg4SEYlhKgIiIjFMRUBEJIapCIiIxDAVARGRGKYiICISw1QERFpoz8Fydh885nUMkZBQERBpgT0Hyzl3ztt8ec4y9h85zp9X5lFRVeN1LJFWs2haVCYtLc3pZjHxUkpG1hfaunXqwJxrTqd/j0QmDuvtQSqR5pnZaudcmr9tOhIQCVDDL0yXn9bvxO9Hjldx14K1fOupldz/t3VU10TPFysRFQGRAG0oPATAzV8eyh9v/hIrfnYp900b0+g5f/loJyPuW8KfV+ZxrKLai5giLaLhIJEA/c9bW3ns7W2suv9y+nbrdKJ976Fyjh6v4ujxar72xPuN9pl16Uh+/NXRxMdZuOOKnKDhIJEQeHdLMROSezUqAAD9eyQyPKkbpw/uyboHJzfa9sQ7OYy4bwlTH1uhYSKJSCoCIgGornFs2n2ICUNOafZ53RM7kpeZTu6j0/j+RcNPtG/afYgR9y0h87XNKgYSUVQERALw9qa9HK+q4azkXgE9Py7OmD3tNDY/NIWffHX0ifa5721nxH1LOOehNyk5WtFGaUUCpyIgEoA1Ow/QMd6YMn5Ai/ZL7BjPXV8Zxc+/NpZR/bqdaN9/tIKzH3qTmc9n8/n+o6GOKxIwTxeVMbM44CGgB5DtnHvOyzwiTcktPsLQPl3pGN+67023nD+MW84fRkVVDTOeW8WKbfsAWLpxL0s37uW756Vw37TTSOig72USXq3+xJnZs2ZWZGbrfdqnmNkWM8sxs4yTvMxVwCCgEihobRaRtpa3/yjD+nYN+nUSOsTx5xmTyH10Gv87Y9KJ9vkf5jH6gdf4xh8+5FB5ZdDvIxKoYL52zAemNGwws3jgSWAqMBa43szGmtnpZvYPn59+QCqw0jl3D3BHEFlE2kxNjWNnSRlDe3cJ2WvGxRkXjOrLw1ePb9S++vNSznhwKSkZWTpnIGHR6uEg59xyM0vxaZ4I5DjncgHM7EXgKufcHGC672uYWQFQ/0n3e2eNmc0EZgIMGTKktXFFWq3o8HHKK2sYGoIjAV/fOXco3zl3KJXVNVz5xAds2n3oxLazH3oTgM0PTSGxY3zI31sEQn9ieBCQ3+BxQV1bU14GrjCz3wHL/T3BOTfPOZfmnEtLSkoKXVKRAOXVnbhN6RO6IwFfHePjeO2HF7JjzjTuumxko21j/uN1Uh94jX1HjhNNN3dKdAh1EfB3W2STn1rnXJlzboZz7i7n3JMhziISErnF9UUg9EcCvsyMn0xOZd2Dk0ns+K//nserakh7+C2+91w2B8o0TCShE+oiUAAkN3g8GCgM8XuIhNX6woN0T+zA4FM6h+09uyd2ZPNDU/kg47JGk9W9vbmIs375JikZWRw5XhW2PNJ+hboIrAJGmdkwM0sArgMWh/g9RMJqw66DjB/YE7Pwz/8zqFdn/njzl8h5ZOoXto3/+Rt8a+5KrWcgQQnmEtEFwEog1cwKzGyGc64KmAW8AWwCFjrnNoQmqkj4Vdc4Nu05zPhBPTzN0SE+jh1zpvHOTy/hv75xxon2j/NKGP3Aa3zpkbc0HYW0SjBXB13fRPsSYEmrE4lEkKLD5VRU1TA0DOcDTsbMGNa3K8P6duXysf1PXD0EUHz4OCPuq/1vt/XhqbrpTAKmT4pIMwoPlAO1wzKRpHfXBPIy01l+76VMGNKr0bbRD7zG7Jc/49P8A55kk+ji6bQRIpEuv6QMgIERVgTqDenThb/deT5Fh8r5/v+uZu3OAwAs+DifBR/XXq2t+wykOToSEGnG2p2ldEmIZ0SS98NBzenXI5G/3Xk+2x+dxtk+RwZj/uN1UjKyKNUdyOKHioBIM3KKj5A6oDsdWjlxXLjFxxkv33k+f7vzvC9sm/BQ7aWlK7YVe5BMIlV0fLJFPFJ4oDxih4KaM2HIKeRlprPuwcmkDW28EM5Nz3xMSkYWf83Ob2JviSUqAiJNcM5ReOAYA3smeh2l1bondmTRHefxQcZl/GxKaqNt9y76jJSMLM1aGuNUBESasOvAMY5X1ZDSBhPHhdugXp2585KR5GWm84srxzXaVj9r6cFjKgaxSEVApAlb9x4GILV/d4+ThNbN56Ww6ZdT+N4Fwxq1n/mL2mKQt08rncUSFQGRJmzdewSAUe2sCAB0Tojngelj2eZnOopLfv0uKRlZ5BYf8SCZhJuKgEgTPsjZR/8enejZuaPXUdpMx/g48jLTyX102he2Xfab90jJyOIfnxVqCut2TEVAxI8DZRWs2LaPK8a1bGH5aBUXZ+yYM41V91/OFeP6N9o264W1DJu9hH1HjnuUTtqS7hgW8aOg9BgA543o63GS8DEzkrp34qmb0gBYuX0/1z/9zxPb0x5+C4C//+B8zkzu5UVEaQM6EhDxo74IRNqcQeH05RF92P7oNP56+5cbtV/15AekZGRRXul3RViJMioCIn4UH66dOK5/z04eJ/FWfJzxpZTe5GWmM+ea0xttq5+Oov4qKolOKgIifpSW1V4zf0qXBI+TRI7rJw4hLzOdH1w6olH75P9eTkpGFnsPlXuUTIKhIiDiR2lZBd07daBjlMwZFE73XjGGvMx07rpsZKP2SY++rTuQo5Cnn3AzG2Jmi83sWTPL8DKLSL3qGsefPsjjsNbwbdZPJqeSl5nOfdPGNGqvvwP5g5x9HiWTlghmeclnzazIzNb7tE8xsy1mlhPAH/bRQJZz7lZgbGuziITSpt2HvI4QVWZeNIK8zHSm+FxOe+MfPyIlI4uNherPSBbMJaLzgSeA5+sbzCweeBL4KlAArDKzxUA8MMdn/1uBtcD9ZvZt4M9BZBEJmfo/WrOnjjnJM6WhuTedA8Df1hbw45c+PdE+7fEVAGTdfQHjBvb0JJs0rdVHAs655UCJT/NEIMc5l+ucqwBeBK5yzq1zzk33+SkCbgF+7py7DEj39z5mNtPMss0su7hY86BL29tQeJCuCfHcduFwr6NEpa9PGExeZjq3X9z4BHL64++TkpHFmxv3epRM/An1OYFBQMNJygvq2pryOnC3mc0F8vw9wTk3zzmX5pxLS0pKCllQkaZs3H2I007tQVyceR0lqmVMrT2B/H93NL7P4Lbns0nJyCLrs90eJZOGQl0E/P2vaXLSEefceufctc65251zPw1xFpEWq6lxbCw8xLiBPbyO0m6cM7T2PoN5dcNF9X7wwhpSMrJ0DsZjoS4CBUByg8eDgcIQv4dIm/m8pIyjFdUau24Dk8cNIC8znbnfaVwMpj62QjedeSjURWAVMMrMhplZAnAdsDjE7yHSZupPCo/VkUCbmTK+thhk3X1Bo/b6m852HTjmUbLYFMwloguAlUCqmRWY2QznXBUwC3gD2AQsdM5tCE1Ukba3ofAgHeKMUf27eR2l3Rs3sCd5melfmJvo/MxlpGRk8aHuMwgLi6Z5wtPS0lx2drbXMaQdu/nZj9l7qJzXf3SR11FiTkFpGRf86p0vtL91z8WM7KeiHAwzW+2cS/O3TffEizSwcfchnQ/wyOBTupCXmc6LM89t1H75b2sXt3ltna4magsqAiJ1ig6XU3z4uK4M8ti5w/uQl5nOkrsvbNR+x19qryZ6fmWeN8HaKRUBkTobdFI4oowd2IO8zHTeuufiRu3/+fcNpGRk8ZZuOgsJFQGROroyKDKN7NeNvMx0lv2kcTH4Xt1NZ9l5vhMXSEuoCIjU2Vh4iCG9u9Ajsf0uLB/NhifVFoN3fnpJo/Zr564kJSOLP67I9SZYlFMREKmzofAgY0/VUUCkG9a3K3mZ6Wz4xRWN2h/O2kRKRhbLNmuYqCVUBESAI8eryNtfppPCUaRrpw7kZabz6c8nN2q/dX7tMNHfP9nlUbLooiIgwr/WEBg3SEUg2vTs3JG8zHQ2/rLxkcEPX/xEs5YGQEVABNiw6yAAY0/VPQLRqktC7ZHB0h83vtGvftbSV9bqyMAfFQERam8S69M1gf49OnkdRYI0un938jLT+fi+rzRq/9FLtUcGf83Ob2LP2KQiIELtPQJjB/bATGsItBf9eiT6PWdw76LPSMnIYvGnmuAYVAREqKlxbCs6wpgB3b2OIm2g/pyB79VEdy9YqyMDVARE2Hu4nIqqGob26ep1FGlD9VcT+Z5Arj8yiNX7DFQEJObt2HcUgKF9unicRMKh/gTyloenNGqvv8/gt0u3eJTMGyoCEvMWfJxPQoc4UjUcFFM6dYgnLzOddQ82Pmfw+LIcUjKy+O83t3qULLzCVgTMbLiZPWNmixq0dTWz58zsaTO7MVxZROo553h/WzFXnjmQft0TvY4jHuie6P+cwWNvbyMlI4v7/7bOo2ThEVARMLNnzazIzNb7tE8xsy1mlmNmGc29hnMu1zk3w6f5GmCRc+424MoWJRcJgYLSY5SWVXJWci+vo4jH6s8ZbH90WqP2v3y0k5SMLDJf2+xRsrYV6JHAfKDRAJqZxQNPAlOBscD1ZjbWzE43s3/4/PRr4nUHA/Wn5qtbHl8kOLl15wNG99dQkNSKjzPyMtPZMadxMZj73vZ2eWQQUBFwzi0HfOdrnQjk1H3DrwBeBK5yzq1zzk33+Slq4qULqC0EAWcRCaX8kjIAhvTWSWFpzMz8njOoPzJ46r3tHiULrWD+8A7iX9/iofYP+qCmnmxmfcxsLjDBzGbXNb8MfMPM/gC82sR+M80s28yyi4uLg4gr8kX5JWUkdIijX3fdKSz+1Z8zyHlkaqP2Oa9tJiUji7sWrPUoWWh0CGJff7dWNrlqvXNuP3C7T9tR4Jbm3sQ5Nw+YB7ULzbc8pkjTthcfIaVPF+LidKewNK9DfBx5mek45xg2e8mJ9lc/LeTVTwv5/sXDmT31NA8Ttk4wRwIFQHKDx4MB3YctUWXznsOkDtDMoRK4+mGiXJ8TyE+9l0tKRhb3LPzEm2CtFEwRWAWMMrNhZpYAXAcsDk0skbZXfPg4BaXHOO1UnRSWlourO4Gcl5neqP3lNbtIycjipmc+8ihZywR6iegCYCWQamYFZjbDOVcFzALeADYBC51zG9ouqkhoLd24B4BLU5u6eE0kMP6uJlqxbR8pGVl854+RXQzMuegZZk9LS3PZ2dlex5B24qZnPiK/pIx3fnqJZg+VkKmpcQy/b8kX2qedPoDf33iOB4nAzFY759L8bdNlmRKTqqpr+HhHCZeN6a8CICHV1DDRknV7SMnI4iu/edebYE1QEZCYlLf/KMerahiv5SSlDfkrBtuLj5KSkcW1f/jQo1SNqQhITPp4RymAJo2TsPBXDLI/LyUlI4vTH3zDo1S1VAQkJj2/Mo+xp/bgNF0eKmHkrxgcLq/ydJhIRUBiTmV1DduLj3BxapJuEhNPNDdMdPZDb4Y1i4qAxJz8kjIqqx0jkrp5HUVinL9iUHK0gpSMLIbNzgpLBhUBiTnbi2tnDh2RpOUkJTL4KwbOQUpGFikZbVsMVAQk5uQWHwFguI4EJML4KwZQWwweydrYJu+pIiAxZ3vxEfp260TPzh29jiLil79i8PSKHW3yXioCEnM27znMyH4aCpLI17AYtNU9LSoCElMOllWyftdBJg7r43UUkYCNG9iDAT3aZg1sFQGJKR9u30eNgwtH9fU6ikhEUBGQmLIiZx/dOnXQwvIidVQEJKa8v20f5w7vTcd4ffRFQEVAYsjO/WXsLCnjgpEaChKppyIgMeOtTXsBuESLyIicENYiYGbDzewZM1vUoO1qM3vazP5uZpPDmUdiy9KNe0jt352Uvro8VKRewEXAzJ41syIzW+/TPsXMtphZjpllNPcazrlc59wMn7ZXnHO3Ad8Fvt2C7CIBKzlawcc7Spg8rr/XUUQiSocWPHc+8ATwfH2DmcUDTwJfBQqAVWa2GIgH5vjsf6tzrqiZ13+g7rVEQi47r4QaB5ekJnkdRSSiBFwEnHPLzSzFp3kikOOcywUwsxeBq5xzc4Dpgbyu1a7tlwm85pxb42f7TGAmwJAhQwKNK9LIhsJDxBmMPbWn11FEIkqw5wQGAfkNHhfUtfllZn3MbC4wwcxm1zXfBVwOXGtmt/vu45yb55xLc86lJSXpW5y0zobCgwzr25XOCfFeRxGJKC0ZDvLH34ocrqknO+f2A7f7tD0OPB5kDpEmlVdW8+H2/Vw9ocnvJyIxK9gjgQIgucHjwUBhkK8pElIrc/dTVlHNlHEDvI4iEnGCLQKrgFFmNszMEoDrgMXBxxIJnTWflxIfZ6SlnOJ1FJGI05JLRBcAK4FUMyswsxnOuSpgFvAGsAlY6Jzb0DZRRVpnzc5SxgzoTpeEYEc/RdqfllwddH0T7UuAJSFLJBJCzjk2FB5i6ngNBYn4o69G0q4VHz7OgbJKUvt39zqKSKv917VnkNBGkx6qCEi7tmH3IQBGD1ARkOg1bmDb3d+iCeSkXfvt0q307prA6YN0k5iIPyoC0m4dPV7Ful0HufnLKXRP1KLyIv6oCEi7tWPfUQBG9e/mcRKRyKUiIO3Wu1tq5ysc2U9FQKQpKgLSLu3Yd5RfL91Kl4R4hvbp4nUckYilIiDt0p8+2EFCfByv/OB8OnXQpHEiTVERkHZpxbZ9XDS6L6N1f4BIs1QEpN3Zc7CcHfuOcu7wPl5HEYl4KgLS7ixZtxuAC0dp/QmRk1ERkHZn0eoCzhjck1TdJSxyUioC0q7sPniMjbsPMf2MU72OIhIVVASkXVm2ufbegEtT+3mcRCQ6qAhIu/LK2l2MSOqqG8REAhS2ImBmw83sGTNb5NPe1cxWm9n0cGWR9mnPwXJW5ZVyzdmDMfO3/LWI+AqoCJjZs2ZWZGbrfdqnmNkWM8sxs4zmXsM5l+ucm+Fn078DCwOPLOLfmp2lAJw/sq/HSUSiR6DrCcwHngCer28ws3jgSeCr1C44v8rMFgPxwByf/W91zhX5vqiZXQ5sBBJbnFykgZoaxytrd5EQH8dpp+qqIJFABVQEnHPLzSzFp3kikOOcywUwsxeBq5xzc4BAh3YuBboCY4FjZrbEOVcT4L4iJyxYtZOlG/dyztBTNE2ESAsEc05gEJDf4HFBXZtfZtbHzOYCE8xsNoBz7n7n3I+AF4Cn/RUAM5tpZtlmll1cXBxEXGnPsj7bTUKHOH53/QSvo4hElWCWl/R35s019WTn3H7g9ia2zW9mv3nAPIC0tLQmX19iV+nRCj7aUcLtFw9nYK/OXscRiSrBHAkUAMkNHg8GCoOLI9JyyzYXUV3jmDx2gNdRRKJOMEVgFTDKzIaZWQJwHbA4NLFEArd04x4G9EjUOsIirRDoJaILgJVAqpkVmNkM51wVMAt4A9gELHTObWi7qCJfdKyimve2FjN5XH/i4nRvgEhLBXp10PVNtC8BloQ0kUgLvLuliPLKGg0FibSSpo2QqFV6tIL/+Pt6knt3ZtLw3l7HEYlKwVwdJOKp7M9L2Xekghdum0THeH2fEWkN/c+RqLV172EAnRAWCYKOBCQq/dfrm/n9u9sZ2DOR7okdvY4jErV0JCBR6ffvbgfgwLFKj5OIRDcVAYk6VdX/ml1k5kXDPUwiEv00HCRRJ7/0GACZ15zOdROHeJxGJLrpSECiTk7REQBGayF5kaCpCEhUqaiq4cWPdxJnaAlJkRBQEZCoUVPjuPMva3h7cxH/OX0sPXRVkEjQVAQkamwoPMRbm/Zy7xWpfPf8YV7HEWkXVAQkaqzIqV1U6Jtpgz1OItJ+qAhI1FixdR9jBnSnX3ctSS0SKioCEhXmvredlbn7uWKcZgsVCSUVAYl4NTWOectzuXh0End/ZZTXcUTaFRUBiXgbCg9RcrSCqycMJF4Lx4iEVNjuGDaz4cD9QE/n3LV1bXHAQ0APINs591y48kj0WL6t9oTwhaOSPE4i0v4Eurzks2ZWZGbrfdqnmNkWM8sxs4zmXsM5l+ucm+HTfBUwCKikduF6kUaqaxz/t6aAs5J70bdbJ6/jiLQ7gQ4HzQemNGwws3jgSWAqMBa43szGmtnpZvYPn59+TbxuKrDSOXcPcEfr/gnSXlVU1fDkOznkFh/ltgs1UZxIWwh0jeHlZpbi0zwRyHHO5QKY2YvAVc65OcD0AN+/AKio+73a3xPMbCYwE2DIEE0WFkt+s3QLTy3PZWJKb6aM11VBIm0hmBPDg4D8Bo8L6tr8MrM+ZjYXmGBms+uaXwauMLPfAcv97eecm+ecS3POpSUlaUw4VpRXVvNSdj5Txg3gpe+fqxPCIm0kmBPD/v5Xuqae7JzbD9zu01YG+J4nEOGNDXs4UFbJDZOGYKYCINJWgjkSKACSGzweDBQGF0ek1p8+yGNony5cMLKv11FE2rVgisAqYJSZDTOzBOA6YHFoYkksW7OzlE/yD3DLeSnEaRhIpE0FeonoAmAlkGpmBWY2wzlXBcwC3gA2AQudcxvaLqq0d/uPHOej3P088/4Ouid24JtpySffSUSCEujVQdc30b4EWBLSRBKzHlmyiZfX7CLO4LYLh9O1k1Y/FWlrmjZCIkJNjWP51n2YQXyc8W/npXgdSSQm6KuWRITNew6z78hxHvn6eC4Y2ZdBvTp7HUkkJqgISESonx/o8tP607+H1gsQCRcNB0lEeGvjXsYM6K4CIBJmKgLiqUPllRSUlpH9eSnTzzjV6zgiMUfDQeKZgtIypj624sSUEF87c6DHiURij4qAeMI5x/1/W09ldQ3HK+Gs5F4M7dPV61giMUdFQDzx4fb9vLe1mP+cPpaLU5Po3DHe60giMUlFQDzx0qp8eiR24IZJQ0hUARDxjE4MS9gdLKvk9Q17uHrCIBUAEY+pCEhYHauo5pElG6moquFbmhtIxHMaDpKwWFdwkHsWfkLe/qNUVjtuPX8Y4wb28DqWSMxTEZA2t23vYf7t2Y/oktCBGRcM59LUJCYN7+N1LBFBRUDa2M79Zdz4x4/oGB/HC7dN0mWgIhFGRUDaxKHySl5bt5vfLcuhorqGhd//sgqASARSEZCQO3K8iit/9z55+8sY1rcrz90ykdH9u3sdS0T8CFsRMLPhwP1AT+fctXVtQ4AngH3AVudcZrjySNv55asb2FlSxrPfTePS1H5aKF4kggW6vOSzZlZkZut92qeY2RYzyzGzjOZewzmX65yb4dM8Gshyzt0KjG1Rcoko+SVl/OCFNdzw9D9ZmF3AHZeM4LIx/VUARCJcoPcJzAemNGwws3jgSWAqtX/ArzezsWZ2upn9w+enXxOvuxa4zsyWAe+07p8gXquqruGuBWt5Z3MRZRXVfCttMD/8ymivY4lIAAJdY3i5maX4NE8EcpxzuQBm9iJwlXNuDjA9wPe/Bfh53esvAv7k+wQzmwnMBBgyZEiALyvh9MQ7OXySf4AnbpjA9DM0E6hINAnmnMAgIL/B4wJgUlNPNrM+wCPABDObXVcsXgceNLMbgDx/+znn5gHzANLS0lwQeSWEVuWVcPufV1NSVoFzcM2EQSoAIlEomCLgb7C3yT/Szrn9wO0+beuBa4PI4Kl3txQx973t3HL+MCaPjZ3x7/1HjnPXC2vp0imeGyaNpEdiR26YpKM0kWgUTBEoABpO/jIYKAwuTvQoKC3j7gVrKauo5p+5JZw9pBf/PmVMu78TtqbGcc/CTykpq+DlO85j/KCeXkcSkSAEUwRWAaPMbBiwC7gOuCEkqSJcZd2JUOdg6Y8v4uMdJfzPW9v49rx/cmlqEj+bMobTTm0f8+IcKq/k0/wDrN15gDU7S/kk/wAHyip5+OrxKgAi7UBARcDMFgCXAH3NrIDak7nPmNks4A0gHnjWObehzZJGkF+/sYW1Ow/w5A1nMzypG8OTunH1hEHM/zCP37+Tw7THV3D1WYO456ujSe7dxeu4Lfb5/qMsWbeHZZv3svrzUmocmMGoft24YuwALhjVV+sBi7QT5lz0nGtNS0tz2dnZnmZYtnkvt87P5jvnDuHhq0//wvaDZZX84b3t/OmDHdQ4x42ThjLrspH07dbJg7Qtt2JbMd97LpvjVTWMG9iDy8b0Y+Kw3pyZ3IseiR29jicirWBmq51zaX63qQgEbvfBY0x7bAWn9uzMy3ee1+yCKHsOlvPY21tZmF1AYoc4vnfhcG67aDjdOjV/8FVRVcPug8fonBBPv+6Jof4nALXX9S9Ylc/iT3ZxxbgB3DhpKJ0T4nl/2z5mPLeKYX278seb0xh8SvQdxYjIF6kIhEBVdQ3XP/1PNhYe4tW7LmB4UreA9ttefITfLN3CknV76NM1gVmXjWTyuAEUHjhGfkkZ+SXHyC8tY2dJGQUlZew5VH5i+GViSm+mnzmQaeMH0CcERxLOOd7dWswjWZvIKTrCoF6d2XXgGH27JfCNswcz/8M8hvXtygu3nUvvrglBv5+IRAYVgRD49RtbeOKdHB677iyuOmtQi/f/JP8Av3ptMytz939hW/8enUg+pQvJvbuQfEpnBvfuwu4D5Sz+dBfbi48SH2ecN6IPXztzIFeMG0DPzi0flik9WsFP/vopyzYXkdKnC7Onncbksf3J/ryUx97axvs5+xgzoDt/+d6kkBQcEYkcKgJBWrGtmH979mO+nZZM5jfOaPXrOOdYuX0/2/cdJfmUziT37sKgXp2bHFZyzrF5z2Fe/bSQVz8rJL/kGAnxcUwe158fXT6Kkf0Cm5lz855D3PZ8NnsPHufeK1K5+bwUEjo0njFkY+EhBvfurHF/kXZIRSAIRYfKmfrYCvp268QrPzifzgneLIzunOOzgoP8/ZNCXlq1k2OV1Vw9YRA/+spohvRpeuz+9fW7uWfhp3Tr1IGnbjqHCUNOCWNqEYkEzRUBrSfQjOoaxw9f/ISyimqeuGGCZwUAwMw4M7kXZyb34geXjmDue9t5fuXnLP6kkG99KZm7LxvFgJ6JlFdWs/tgObsPHOO9bcU89V4uZyX34qmbzqF/j7Y50Swi0UtFoBm/W7aNlbn7+X/XnsGoCFoUpU+3TtyfPpbvXTic3y3bxkur8lm0uoCuCfGUllU2eu615wzm4avHN3slk4jELhWBJny4fR+Pvb2Na84exDfTkk++gwf690jk4atP5/sXjeCZ93dQUV3DwJ6JnNqzM6f2SmRwry7NDhWJiKgI+LHvyHF+9OInDO/blYeuGu91nJNK7t2FB68c53UMEYlCKgI+amocP37pEw4eq+T5GRPpepKbu0REopn+wvn4w3vbWbFtH3OuOZ0xA9rHJHAiIk0JdHnJmLAqr4TfvrmVr505kOu+FJnnAUREQklFoE7J0QruemEtyad05tGvj4+ZBWJEJLZpOIja8wA//eunlByt4OU7z6O77poVkRihIwHgmfd3sGxzEQ9MP00LpYhITAnbkYCZXQ2kA/2AJ51zS82sK/B7oAJ41zn3l3DlqbdmZym/en0zU8cP4KZzh4b77UVEPBXQkYCZPWtmRWa23qd9ipltMbMcM8to7jWcc684524Dvgt8u675GmBRXfuVLY8fnINlldz1wlpO7ZVI5jfO0HkAEYk5gR4JzAeeAJ6vbzCzeOBJ4KvULjq/yswWU7vU5Byf/W91zhXV/f5A3X5Quzj9urrfq1saPhjOOe5d9ClFh8tZdPt5rZqeWUQk2gVUBJxzy80sxad5IpDjnMsFMLMXgaucc3OA6b6vYbVfszOB15xza+qaC6gtBJ8Q5vMT8z/MY+nGvTyQfhpnJvcK51uLiESMYP7wDgLyGzwuqGtryl3A5cC1ZnZ7XdvLwDfM7A/Aq/52MrOZZpZtZtnFxcVBxP2XzwoO8OiSTVx+Wn9mXDAsJK8pIhKNgjkx7G8AvcnFCZxzjwOP+7QdBW5p7k2cc/OAeVC7nkDLYzZ2qLySWS+sJalbJ379TZ0HEJHYFkwRKAAa3lY7GCgMLk7bcs4x++V17DpwjIXfP5deXbSOrojEtmCGg1YBo8xsmJklANcBi0MTq2385aOdZH22m3uvSOWcob29jiMi4rlALxFdAKwEUs2swMxmOOeqgFnAG8AmYKFzbkPbRQ3OhsKD/PIfG7kkNYmZFw73Oo6ISESIqjWG+48Y62589IVW7fvRjhKqampYcveF9OnWKcTJREQiV7tZY7jseDWrPi9p1b7dOnXgka+PVwEQEWkgqopA6oDurPjZZV7HEBFpNzSBnIhIDFMREBGJYSoCIiIxTEVARCSGqQiIiMQwFQERkRimIiAiEsNUBEREYlhUTRthZsXA5w2aegIHfZ7WsM13e19gX4hj+csQin2ae05T207WH76Pw9E/TeUK9vnh6B/fx9HUPyd7XiB94a8tGvon0H1iqX+GOueS/G5xzkXtDzCvuTbf7UB2ODKEYp/mntPUtpP1hxf905o+ipT+8dNfUdM/remj9tI/bfUZak/90/An2oeD/K1G9upJtocjQyj2ae45TW07WX/4Pg5H/7TmfSKlfwLNEqy26J+TPS+QvvDXFg39E+g+sdw/J0TVcFCwzCzbNTGTnqh/Tkb90zz1T/MitX+i/UigpeZ5HSDCqX+ap/5pnvqneRHZPzF1JCAiIo3F2pGAiIg0oCIgIhLDVARERGJYzBYBM+tqZs+Z2dNmdqPXeSKNmQ03s2fMbJHXWSKVmV1d9/n5u5lN9jpPpDGz08xsrpktMrM7vM4Tier+Dq02s+leZWhXRcDMnjWzIjNb79M+xcy2mFmOmWXUNV8DLHLO3QZcGfawHmhJ/zjncp1zM7xJ6p0W9tErdZ+f7wLf9iBu2LWwfzY5524HvgVE3KWRbaGFf4MA/h1YGN6UjbWrIgDMB6Y0bDCzeOBJYCowFrjezMYCg4H8uqdVhzGjl+YTeP/Eqvm0vI8eqNseC+bTgv4xsyuB94G3wxvTM/MJsH/M7HJgI7A33CEbaldFwDm3HCjxaZ4I5NR9s60AXgSuAgqoLQTQzvqhKS3sn5jUkj6yWr8CXnPOrQl3Vi+09DPknFvsnDsPiIkh1xb2z6XAucANwG1m5snfoQ5evGmYDeJf3/ih9o//JOBx4AkzSyd80ydEIr/9Y2Z9gEeACWY22zk3x5N0kaGpz9BdwOVATzMb6Zyb60W4CNDUZ+gSaoddOwFLwh8rYvjtH+fcLAAz+y6wzzlX40G2mCgC5qfNOeeOAreEO0wEaqp/9gO3hztMhGqqjx6n9stErGuqf94F3g1vlIjkt39O/OLc/PBF+aJYGAYpAJIbPB4MFHqUJRKpf05OfdQ89U/zIrp/YqEIrAJGmdkwM0sArgMWe5wpkqh/Tk591Dz1T/Miun/aVREwswXASiDVzArMbIZzrgqYBbwBbAIWOuc2eJnTK+qfk1MfNU/907xo7B9NICciEsPa1ZGAiIi0jIqAiEgMUxEQEYlhKgIiIjFMRUBEJIapCIiIxDAVARGRGKYiICISw1QERERi2P8HXhypYODDoDIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(hp.anafast(map_log_pol_tens_small_scale[0], lmax=16384))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combine scales\n",
    "\n",
    "* Combine small and large scale maps\n",
    "* Transform from logpoltens to IQU\n",
    "* Write output map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_log_pol_tens = map_log_pol_tens_large_scale\n",
    "map_log_pol_tens += map_log_pol_tens_small_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "del map_log_pol_tens_small_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pysm3.utils import log_pol_tens_to_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_map = log_pol_tens_to_map(map_log_pol_tens)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.72037508e+01,  3.66491306e+01,  3.67077472e+01, ...,\n",
       "         4.10425484e+01,  4.07681616e+01,  4.11603859e+01],\n",
       "       [ 7.03200141e-01, -6.66901679e-01,  6.07019814e-01, ...,\n",
       "         7.62817241e-01, -6.98765800e-01,  7.95575034e-01],\n",
       "       [-8.09556476e-02,  1.21250909e-01,  8.08351243e-03, ...,\n",
       "        -8.87983650e-01,  8.63559839e-01, -9.42672655e-01]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "output_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7faae20f1460>]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAl3ElEQVR4nO3deXhV1bnH8e+bGcIMYQpDwiCjoBIRUREUNRgc6lBFW6+WSrVqW3tvNShaWkVT7YBWq0WrqG1BqtYJUARBBkEIggwyRQgQpgQiMwSSrPtHYsxMQnLOPsn5fZ6HR/Z7hv2e9cRfNuvsvbY55xARkfovxOsGRETEPxT4IiJBQoEvIhIkFPgiIkFCgS8iEiQU+CIiQSLM6wYq06pVKxcXF+d1GyIidcry5cv3OudiStcDOvDj4uJITU31ug0RkTrFzLaWV9eUjohIkFDgi4gECQW+iEiQUOCLiAQJBb6ISJDw21k6ZhYCPAY0AVKdc6/5a98iIlLDI3wze8XMMs1sTal6opltMLM0M0suLF8DxAIngYya7FdERKqvplM6k4HE4gUzCwWeB0YAvYFRZtYb6AEsds79Gri7hvsVEZFqqlHgO+fmA9mlygOBNOfcZufcCWAqBUf3GcC3hc/Jq+g9zWyMmaWaWWpWVlZN2hMRkWJ88aVtLLC92HZGYe0d4Aoz+yswv6IXO+cmOecSnHMJMTFlrgwWEZHT5Isvba2cmnPOHQVG+2B/IiJSBb44ws8AOhbb7gDs9MF+RESkGnwR+MuA7mYWb2YRwM3A+z7Yj4iIVENNT8ucAiwGephZhpmNds7lAvcCHwPrgGnOubU1b1VERGqiRnP4zrlRFdRnADNq8t4iIlK7tLSCiEiQUOCLiAQJBb6ISJBQ4IuIBAkFvohIkFDgi4gECQV+NaRlHubA0ZNetyEiclr8dgOUuuyr7ft5Zs4mPl2fyZmxTXnn54MJD9XvShGpW5RalVix7Vtuf3Up1zy/iC+3fcuNAzqwescBJs7e6HVrIiLVpiP8cizfms3E2ZtYsGkvzRuG80BiD247P45GkWGYwd/mfcPFZ7RmYHwLr1sVEakyBX4xy9KzeWb2Jham7aVldATJI3ry40GdiY78fpgevaoPX2zJ5v43VzLzVxfRJCrcw45FRKpOgQ8s2byPZ2ZvYvHmfbRqFMHDV/bi1kGdaBhRdngaRYbxl5vO4sYXFzP+vbX8+aaz/N+wiMhpCNrAd86xuDDov9iSTUzjSMYl9eLW8zrTICK00tee06k59w7rxjNzNjGsZ2uu6t/eT12LiJy+oAt85xyff1MQ9EvTs2ndOJJHR/bmlvM6ERVeedAXd98l3fhsYxYP/3c1Azo3p32zBqfdU36+IySkvBuFiYjUnqA5S8c5x/yNWdz44mJuffkLtmUf5XdX92H+A8P4yYXx1Qp7gLDQECbedBa5+Y7/nfYV+fnutPpauiWbAY9/wqy1u0/r9SIiVVXvA985x7wNmVz3wufc9spSduw/xmPX9GHeb4byP4Pjqh30xcW1imb8VX1YvHkfLy/cXO3Xb9t3lJ+9kcq3R0/yx1kbTvuXhohIVdTbKZ2CoM/imTmbWLl9P7HNGvD4tX25MaEDkWGnH/Kl3ZjQgTnr9/D0xxu4oFsr+rRvWqXXHTp+ktGvLSPfwa+Gd2fi7E18sm4PV/RpW2u9iYgUZ84F7lFlhzP6ul88+9ZpvXZh2l5WZRygQ/MG3DOsG9ef04GIMN/8gyb7yAkSJ86naYNwPrjvwlP+qyEv3/HT15axYNNeXh89kIFxLbjkT5/RrGE4791zAWaazxeR02dmy51zCaXrAX2Ev/dQDi8tqP5UCUBsswb84fozue6cDj5fBqFFdAR/vLE/t72ylJSZ6xl/dZ9Kn//kjHXM3ZDFhB/0ZXDXVgDcPbQrY99ZzYJNexlyRkylr3fOcfREXonrA0RETiWgE6NvbFNSJ1zpdRtVMuSMGO64II5XF6UztEcMQ3u0Lvd5U5du4+WFW7h9cBy3nte5qH7dObE8M3sTz32adsrAf2nBZibO3sTc/xtKmyZRtfo5RKT+qvdf2vrTg4k9OaNNI37z1ir2Hc4p8/iSzfsY9+4ahpwRw7ikXiUeiwwL5WcXd2FpejZLt2RXuI+d+4/xl082cfREHm8tz6j1zyAi9ZcCvxZFhYcy8aazOXD0JGPfWU3x70e27jvC3f9cTlyraJ675WzCyplmuvncTrSMjuC5uWkV7mPCjHXkO0evdk2YumybzuwRkSpT4Ney3u2b8JsrejDr6z28uWw7AAePn2T0a6k44B//k1Dh+jsNIkIZfVE88zdmsSpjf5nHF6XtZfqqXfx8aDfuurgL27OP8fk3+8p9r017DjHu3dWcyM2vrY8mInWcAt8HRl8Yz+CuLfndB1+TlnmYe/+9gvS9R3jh1gF0bhld6Wt/PKgzTaLCeL7UUf7JvHx++/5aOrVoyM8u7sIVfdrSrGE4U5ZtK/d9npy5nn8u2caCTVm19rlEpG5T4PtASIjxpx/2JyIshGueW8j8jVk8fm1fzu/a8pSvbRwVzu2D4/h47R427jlUVH/t83TSMg/z6MjeRIWHEhUeynVnd2DW2t1lvi9Yt+sgn67PBOCDr3bW7ocTkTrLr4FvZtFmttzMRvpzv15o17QBT/zgTI6cyOMnF8Rz88BOVX7tHRfE0zAilL8VHuVnHjzOxNmbGNYjhkt7fX/2z6iBHTmZ53jnyx0lXv/CvG+Ijggl6cx2fPL1Ho6dyCvx+Pbso/xp1gZy8zTdIxJMahT4ZvaKmWWa2ZpS9UQz22BmaWaWXOyhB4FpNdlnXZLUrx0LHxzGIyN7nfrJxTSPjuDW8zrx/lc7Sd97hCdnrudEbj6/vapPiYuyurdpzIDOzZmybFvRF8Rb9x3hw1U7+dGgztw6qBNHTuQxd0Nmiff/8ycb+eunaRXO/4tI/VTTI/zJQGLxgpmFAs8DI4DewCgz621mw4GvgT013Ged0qF5w9O6cvbOi7oQFhrC/dNW8t8VOxgzpAtxrcrO/998bkc2Zx1hWfq3ALz42WbCQkMYfWE858W3JKZxZIlpncxDx/lwVcH2d/8VkeBQo8B3zs0HSp80PhBIc85tds6dAKYC1wDDgEHALcCdZlbuvs1sjJmlmllqVlbwfuHYukkUNyV0ZMW2/bRvGsXPh3Ut93lJ/drRODKMqUu3sefgcd5ensENAzrQukkUoSFG0pnt+HR9JoeOnwTgX0u2kZvvODeuOR+t2V3iLJ4X5n3D7z/42i+fT0T8zxdz+LHA9mLbGUCsc+5h59yvgH8DLznnyp1Ads5Ncs4lOOcSYmIqv+K0vrtraFe6tW7E4z/oW+7dtwAaRoRx9Vntmb56F3/5ZCO5+fncNeT7Xw5X9W9PTm4+s9ftISc3j399sY1hPVpz99CuHDyey8K0gl+q+w7nMHH2RiZ/voW95Vw0JiJ1ny8Cv7z5i6Krg5xzk51zH/pgv/VObLMGzP71xVzSs02lzxs1sBM5uflMXbadq/q3p1PLhkWPndOpGbHNGvDBV7uYsXoXew/ncPvgOC7sFkPTBuF8+NUuAP65ZBs5ufnkO5i5elfR6+9/cyWT5n/jmw8oIn7li8DPADoW2+4AaLLYh/rGNqVvbBOgYBG24syMkf3bMX9jFi/M+4auMdFc1L0VEWEhXNGnDbO+3sOBoyd5fXE6l/RsTffWjfig8JfAqoz9/HfFDl5esEVX9IrUA74I/GVAdzOLN7MI4GbgfR/sR4oZf1Uffnd1H3q2bVLmsav6tSc337Fxz2FuHxxX9CXyyH7tOZyTyy/fXMG+IycYM6QLV/Vvz9L0bHYdOMari9IByDyUw4rtBV8KT5j+NT99bZnfPpeI1J6anpY5BVgM9DCzDDMb7ZzLBe4FPgbWAdOcc2tr3qpUJiGuBf8zOK7cx/q0b0KXVtE0jgzjunM6FNUHd21Ji+gI5m3Iol+HppwX34KR/doBMHlROh+u2skNAzoQERrCzNUFF3i9tngrs9dlsuvAMQDSMg9zOCfX559PRGquRssjO+dGVVCfAcyoyXtL7TEznr6xH0dySq6hHxYawoi+bfnXF9sYM6QLZkaXmEb0jW3C3+dvxgzuHdaN7CMnmLlmN82jI4rO6pm9LpOkM9uR9OwCRg3sdMp7AIiI9wJ6PXypPQM6tyi3/rMhXWkRHUFisVsrjuzXnjU7DnJJj9bEtYpmRN+2fLo+k7/NTeOCbi3Zuf944RW8ueTk5jN3Qybj6cPCTXvJOnycH5zdodx9iYi3FPhBrlPLhvzv5T1K1K49K5Zpqdv5+bBuAFzWuw1hIcaRE3ncPjieZenZvLpoC1v2HiY0xNi67yhb9h5h3LuryTqUw8h+7X1+lzERqT79XylltG0axaf/O5QBnZsD0KxhBBd1b0Vcy4Zc0rM1l/Vuw8k8x/bsY9xT+EvhmdkbSd93lCMn8li94wDPz01j7DurvfwYIlKKjvClSibefDYn8/IJDTHO6dScltERnMzL5+dDu/Leyh28u3InkWEh5OTmM39jFq8uSudwTi7JiT3ZeySHmMaRFd4HQET8Q0f4UiVNG4TTqlEkAKEhxqNX9ebxH5xJVHgoFxfeg3dE37b0bNuYfyzcwoFjJ8nLd3ywaicjn13IUx+t97J9EUGBL6fpmrNiubp/ewCG9yq4EviHCR0Z3LUVh47n0rRBOE0bhJMycz3HTuYxd30W763cwbkTZus0ThGPKPClxoacEcNnvxnK4G6tGFx4k5crz2zLxWfEcDgnl/BQY8f+Yzw5Yz1Zh3JYuCmLX09byaZiN3gREd9T4Eut+O7WjYO7teTy3m24fXA8w3oWTPXcf9kZAOw+eByApz/ewDtf7uA/yzNYlp5NZmFdRHxLX9pKrWoYEcak2xIA6BoTTYPwMC7v3Yb/frmDzXuP0KVVNJsyDwMwa+1uXl20hWvOiuXhK3vRMDKUyLBQL9sXqdcU+OIzYaEhJPYtuKDrF5d2J+PbYxw7kcumT9OICAshfd9RAOZtyGLehiyu6t+OR0f2Pq0bxojIqWlKR/ziqv7tuXtoV4b2bI0Z3Fd4/r4Z7D2cw97DOXy0ZjcDn5jD64vTee3zdI6fzDvFu4pIdegIX/zqnE7NWT7uMqIjQ3l9yVauOzuWv8/fDMCuAwVz+ROmryMnN59jJ/PoF9uUwd1aedmySL2hwBe/axEdAcCSsZcSYnDg2ElaNYrkublpAOQULtCWMrPg3P1JPx5An9imxDZr4E3DIvWEAl88ExpSMFefcn0/nHNs3HOIhhGhvLuy5P1yxryxnPZNo/jrLWdXuAiciJyaAl8Cgpkx6bYEjp/MIyc3n8M5uSzYtLfo8Z0HjnP9C4v5953nMTCuBWFanE2k2sy5wL11XUJCgktNTfW6DfHAsRN5fLYxi2fnbOLrXQdLPDa4a0seurIXfWObetSdSGAzs+XOuYQydQW+BLJjJ/I4ePwkD769inkbsko8Nv6q3txyXmciwnS0L1KcAl/qvPS9R1i/+yB//TSNtTsLjvp7tm3M87eeQ9eYRh53JxI4Kgp8HRpJnRHXKprEvu2Y/ouLeOzavgCs332IERMXMG9DpsfdiQQ+Bb7UST8e1JkvHrqUHw/qzIm8fG5/dRk//scXulhLpBIKfKmz2jSJ4rFr+7LggWGEhhgLNu2l5yMf8eTMdUU3WxeR7ynwpc7r2KIhaRNGcOOAgpun//2zzdzw4uds3XfE485EAosCX+oFM+PpG/uTOm449wzryrpdB7n46Xk8MWOdbrgiUkiBL/VKq0aR/OaKnrx3z4Vc1L0Vk+Zv5uZJi8n49qjXrYl4ToEv9VLv9k14Y/R5PHfL2aRlHubCP8zlyZnrCOTTkEV8zW+Bb2bXmtlLZvaemV3ur/1KcBvZrz3/+ukgoGBuP37sDLZn62hfglONAt/MXjGzTDNbU6qeaGYbzCzNzJIBnHPvOufuBG4HbqrJfkWqY0Dn5kUrcwJc9NRclmze521TIh6o6RH+ZCCxeMHMQoHngRFAb2CUmfUu9pRxhY+L+E3bplF888SVjCi8A9fNk5bQ97cfs+vAMY87E/GfGgW+c24+kF2qPBBIc85tds6dAKYC11iBPwAznXNf1mS/IqfDzHjhRwN46vp+ABzOyeXKZxbw0ZrdumBLgoIv5vBjge3FtjMKa/cBw4EbzOyuil5sZmPMLNXMUrOysip6mshp++G5HVn/WCLPjjqb3HzHXf9cTs9HPiLrUI7XrYn4lC8Cv7w7UDvn3LPOuQHOubuccy9W9GLn3CTnXIJzLiEmJsYH7YlAVHgoV/dvz5xfX1xUO3fCbNIyD3nYlYhv+SLwM4COxbY7ADsreK6Ip1o3iWLThBF0b12w2ubwP89nz8HjHncl4hu+CPxlQHczizezCOBm4H0f7EekVoSHhjDr/iGMGdIFgPOemMOri7Z43JVI7avpaZlTgMVADzPLMLPRzrlc4F7gY2AdMM05t7bmrYr4jpnx0JW9mDqm4Jz9333wNT/8+2JdqCX1im6AIlLK/I1Z3PbK0qLttAkjdA9dqVN0AxSRKhpyRgzTf3Fh0fb5KZ/qtE2pFxT4IuXo074pnydfAkDWoRwe/u+aU7xCJPAp8EUq0L5ZA1Y8chk/GtSJt7/M4P43V2pOX+o0Bb5IJZpHR/BAYk8iw0L474odPDNnk9ctiZw2Bb7IKTSJCmfV+IIFXifO3sQfP97Atc8v0kVaUuco8EWqIDIslMVjL6Fn28Y8NzeNldv388qidK/bEqkWBb5IFbVr2oD37r2ApH7tAPj3F9tYuX2/t02JVIPOwxc5Dduzj3Lpnz7jRF4+fdo34Yo+bfnFpd29bksE0Hn4IrWqY4uGvHL7uQCs3XmQP3+ykWXppVcKFwksCnyR03Rh91Y8dX0/Yps1AODGF7UUgwQ2Bb5IDfzw3I4sSr6kaLXNxd/o1okSuBT4IrVg8k8GAjBhxjry8nWUL4FJgS9SC2KbNWDMkC6s3XmQHuNm8us3V5Kv4JcAo8AXqSX3XdKNIWfEkJvveGfFDn4xdYXXLYmUoMAXqSWNo8KZ9OMBnNWxGQAfrtrFlr1HvG1KpBgFvkgtigoP5d17LuD1wjn9+99c6W1DIsUo8EV8YMgZMVx/TgdWbt9PXPJ03lqewcHjJ71uS4KcAl/ERx5M7FH09//7z1f0Gz9LX+SKpxT4Ij7SukkU6x9LZFxSr6Jal4dm8OBbqzzsSoKZAl/Eh6LCQ/npRV1IHTe8qPZm6naWb/3Ww64kWCnwRfygVaNI1v0+keG9WgNw/Quf8+aybfzo5S94fXE6zjm2Zx/1uEup77RapoifPfbh1/xj4ZYStTsuiOPVRem8+KNzSOzbzqPOpL7QapkiAeKRkb3L1F4tvJnKEzPW+7kbCSYKfBEPbHnySp65+SyiI0JL1LdlHyU3L9+jrqS+05SOiMecc+w+eJyUmet5b+VOBsa3YNrPzve6LanDNKUjEqDMjHZNG/D0Df0BWLolm7jk6Yx9Z5Uu1pJa5bfAN7NoM3vNzF4ys1v9tV+RuiIirOT/jlOWbufa5xd51I3URzUKfDN7xcwyzWxNqXqimW0wszQzSy4sXwe85Zy7E7i6JvsVqa+WPTy8xPbmrCNs3HNIN1aRWlHTI/zJQGLxgpmFAs8DI4DewCgz6w10ALYXPi2vhvsVqZdiGkeSOm44aRNG0KxhOACX/2U+o15agnOOr7bv572VOzzuUuqqGgW+c24+UPrOzQOBNOfcZufcCWAqcA2QQUHoV7pfMxtjZqlmlpqVlVWT9kTqpFaNIgkLDWHFI5dxTqdmRfXfvr+Wa55fxC+nrvSsN6nbfDGHH8v3R/JQEPSxwDvA9Wb2AvBBRS92zk1yziU45xJiYmJ80J5I3WBmvPPzC4q2X1+8tejvL372DZ+u3+NFW1KHhfngPa2cmnPOHQHu8MH+ROq11eMv58zxs0rUUmYWXKB1We82vHRbmbPvRMrliyP8DKBjse0OwE4f7EckKDSOCic9JYm1v7uCkf1KLrvwydd7mLV2NwDLt36r9XikUjW+8MrM4oAPnXN9C7fDgI3ApcAOYBlwi3NubXXfWxdeiZRv4aa9/OgfXxRtP5jYkz98VHDUn56S5FVbEiB8cuGVmU0BFgM9zCzDzEY753KBe4GPgXXAtNMJexGp2IXdW/HuPd/P738X9iKV0dIKInXYjv3HuCDl0zL1mMaRXN67DRN+cKYHXYnXKjrC98WXtiLiJ7HNGhRN4Wzac4jL/jIfgKxDOfzri200CA/l4aRemJV3LoUEG62lI1JPdG/TmAUPDCtRe3nhFuLHziCQ/yUv/qPAF6lHOrZoyIbHE9n8xJU8dGXPonr82BkcO6EL3IOdAl+knokMCyUkxBgzpGuJe+n2evQjPtuoq9eDmb60FanncvPy6fbwzDL1db9PpEGpG7BI/aD18EWC1Hfr8pT2g79p6eVgo8AXCQLNoyNY/1iJhW1Zv/sQccnTycsP3H/lS+1S4IsEiajwUNJTkspcidv1oRm8vGCzvtQNAgp8kSCUnpLE5DvOLdp+fPo6ej36EQeO6paK9ZkCXyRIDe3Rmnd+PrhErf/vZ5F56LhHHYmvKfBFgtg5nZqTnpJUdHctgIET5vDV9v3eNSU+o8AXEVY+ejkXn/H9DYeueX4Rc9dnetiR+ILOwxeRIut3HyRx4oKi7cFdW5IQ14L7h3fXejx1iM7DF5FT6tm2CRsfH1G0/fk3+3h2ziZ2HdC8fn2gwBeREiLCQtjy5JUlaoNTPmXo03N5Y3G6N01JrVDgi0gZZlYm9NP3HeWR99aS9OyCCl4lgU6BLyLlMjPSU5L4PPmSEvW1Ow8SlzxdSy7XQQp8EalU+2YNWPrwpUz/xYUl6vFjZ7B8azZb9h7xqDOpLgW+iJxS68ZR9GnflPSUJIb3alNUv/6FxQz74zwd7dcRCnwRqZaXbhtAXMuGJWrxY2dw5viPWZWxX+EfwBT4IlItZsa83wxj9fjLmfaz84vqh47ncvVzi/jXF9s87E4qo8AXkdPSOCqcgfEt2PLklfwwoUNRfdy7a9h/9ISHnUlFFPgiUiNmxlM39Oc3V/Qoqp31+0+4/C+f8cXmfR52JqUp8EWkVtwzrFuJe+hu3HOYmyYt4YOvdnL8pNbaDwQKfBGpNa0aRbLmd1eUqN03ZQU9H/mIv81L86gr+Y5fA9/MrjWzl8zsPTO73J/7FhH/aBQZRnpKEu/ec0GJ+lMfbeCnr6XqLB4PVTnwzewVM8s0szWl6olmtsHM0swsubL3cM6965y7E7gduOm0OhaROuGsjs1IT0nixR8NKKrNXreH+LEztPSyR6pzhD8ZKHEXZDMLBZ4HRgC9gVFm1tvMzjSzD0v9aV3speMKXyci9Vxi37ZsfqLkujx3TF7GweO6naK/VTnwnXPzgexS5YFAmnNus3PuBDAVuMY5t9o5N7LUn0wr8AdgpnPuy9r7GCISyEJCChZj++ON/Ytq/cbP4uKn52qKx49qOocfC2wvtp1RWKvIfcBw4AYzu6u8J5jZGDNLNbPUrKysGrYnIoHCzLhhQAfevvv7i7W27jtK/NgZbM8+6mFnwaOmgV/eLXAq/HXtnHvWOTfAOXeXc+7FCp4zyTmX4JxLiImJKe8pIlKHDejcgrQJI0rULnpqLmPfWeVRR8GjpoGfAXQstt0B2FnD9xSRei4stOAmKxGh30fQlKXbOeexT3hreYameXykpoG/DOhuZvFmFgHcDLxf87ZEpL4zMzZOGFHiRivZR07wf//5ivixMzzsrP6qzmmZU4DFQA8zyzCz0c65XOBe4GNgHTDNObfWN62KSH303Y1WSotLns7qjAMedFR/WSD/0ykhIcGlpqZ63YaI+MnB4yfpN35WidqCB4bRsUXDCl4h5TGz5c65hNJ1La0gIgGjSVR4maP9i56aS1zydPYdzvGoq/pDgS8iASc9JYkvH7msRG3A47N5c5nW2q8JBb6IBKQW0RElvtAFePDt1byxZCv5+YE7FR3INIcvIgFv3+EcBjw+u0RtUJcWTB1zfgWvCG6awxeROqtlo0g2Pl7yYq0lm7OJS56uo/1qUOCLSJ0QERZCekoSs389pES9y0Mz+HLbtx51Vbco8EWkTunWujHpKUlc3rtNUe26v33O6MnLPOyqblDgi0idNOm2BFYUO5NnzvpM4pKn88aSrR52FdgU+CJSZzWPjmDJ2EtL1B55dw3nPznHo44CmwJfROq0tk2jylystevAceKSp7M567BHXQUmBb6I1AvpKUll7qx1yZ8+49aXl5CTm+dRV4FFgS8i9UZISNmF2Bal7aPHuI/YsPuQR10FDgW+iNQ76SlJrPt9iVtwc8XE+dw3ZYVHHQUGBb6I1EsNIkLLLM3wwVc7iUuezpGcXI+68pYCX0Tqre/W2l//WMmj/T6//Zhpy7ZX8Kr6S4EvIvVeVHgo6SlJDIxrUVR74O1VxCVPJ+tQ8Cy7rMAXkaAx7a7zyyy7fO6E2WQePO5RR/6lwBeRoNIiOoL0lCR6tWtSVBv4xJygWIhNgS8iQWnmLy9i1fjLS9S6PDSDjG+PetSR7ynwRSRoNYkKL7Ps8oV/KLilYn2kwBeRoPbdssuPjOxdol4fp3gU+CIiwOgL48tcpdvloRn16mhfgS8iUkx6ShL3Dz+jRC0ueTrf1IOF2BT4IiKl/HJ4d9b+7ooStUv/9BnX/W2RRx3VDgW+iEg5oiPDSE9JYsyQLkW1L7ftJy55Os7Vzbl9Bb6ISCUeurJXmTV54sfOYM66PR51dPr8GvhmFm1my81spD/3KyJSE9+tydO+aVRRbfRrqXXuC90qBb6ZvWJmmWa2plQ90cw2mFmamSVX4a0eBKadTqMiIl77fOylvH334BK1uOTpfLRml0cdVY9VZS7KzIYAh4HXnXN9C2uhwEbgMiADWAaMAkKBJ0u9xU+AfkArIArY65z78FT7TUhIcKmpqVX+MCIi/jJq0hIWb95Xolb6tE6vmNly51xC6XqVjvCdc/OB7FLlgUCac26zc+4EMBW4xjm32jk3stSfTGAYMAi4BbjTzMrdt5mNMbNUM0vNysqqxkcUEfGfKWMGlXuxViDfWasmc/ixQPEFpTMKa+Vyzj3snPsV8G/gJedcfgXPm+ScS3DOJcTExNSgPRER3xp9YXyZL3SvmDifR95dU8ErvFWTwLdyaqecH3LOTa7KdI6ISF1Q3k1W3liylbjk6ezYf8zDzsqqSeBnAB2LbXcAdtasHRGRuum7m6wUd0HKp6TMXO9RR2XVJPCXAd3NLN7MIoCbgfdrpy0RkbopPSWJKXcOKtp+8bNviEueTm5eubPYflXV0zKnAIuBHmaWYWajnXO5wL3Ax8A6YJpzbq3vWhURqRvO79qSzU+UnNvv9vBMfvaGt2cdVum0TK/otEwRqeuS317F1FI3TPf16Zs1Oi1TREROT8r1/cosxBaXPJ256zP93osCX0TEx75biO2i7q2KandMXub3pRkU+CIifvLG6PP46rcl76Prz7X2FfgiIn7UtEF4mTn8S//0mV+O9hX4IiIeSE9J4q27zi9Ri0uezsHjJ322TwW+iIhHEuJalDna7zd+FiP/usAn+1Pgi4h4LD0licVjLynaXrPjoE8u1FLgi4gEgHZNG/BNsYu1so+eqPV9KPBFRAJEaIjxxxv7A5BzUkf4IiL1WlR4QSwfO5lX6++twBcRCSANwkMBOHZCgS8iUq9FhhUEfk6upnREROq1kMJbS/liYUsFvohIAIkMD6FjiwZEFk7t1KawWn9HERE5bQM6t2DBA5ec+omnQUf4IiJBQoEvIhIkFPgiIkFCgS8iEiQU+CIiQUKBLyISJBT4IiJBQoEvIhIkzBeX79YWM8sCthYrNQUOVHG7FbDXB22V3mdtvaay55T3WFVqwTw+5dVPNWalH/fFGAXy+JRX8/f4VNRXTZ/vj/Epve3l+HR2zsWUqTrn6swfYFJVt4FUf/RQW6+p7DnlPVaVWjCPT1XGo5wxKf38Wh+jQB6fKvzM+Hx8TmeMAmV8yhmvgBif4n/q2pTOB9Xc9kcPtfWayp5T3mNVqQXz+JRXP9WYBfv4lFfz9/iczn4CZXyq2ktNnfY+AnpKpybMLNU5l+B1H4FK43NqGqPKaXwqF4jjU9eO8KtjktcNBDiNz6lpjCqn8alcwI1PvT3CFxGRkurzEb6IiBSjwBcRCRIKfBGRIBE0gW9m0Wb2mpm9ZGa3et1PoDGzLmb2DzN7y+teApGZXVv4s/OemV3udT+Bxsx6mdmLZvaWmd3tdT+BqjCHlpvZSC/2X6cD38xeMbNMM1tTqp5oZhvMLM3MkgvL1wFvOefuBK72e7MeqM74OOc2O+dGe9OpN6o5Pu8W/uzcDtzkQbt+V83xWeecuwv4IRBQpyL6UjUzCOBBYJp/u/xenQ58YDKQWLxgZqHA88AIoDcwysx6Ax2A7YVPy/Njj16aTNXHJxhNpvrjM67w8WAwmWqMj5ldDSwE5vi3TU9NpopjZGbDga+BPf5u8jt1OvCdc/OB7FLlgUBa4RHrCWAqcA2QQUHoQx3/3FVVzfEJOtUZHyvwB2Cmc+5Lf/fqher+/Djn3nfODQaCZsq0mmM0DBgE3ALcaWZ+z6Ewf+/QD2L5/kgeCoL+POBZ4DkzS8J/l4gHonLHx8xaAhOAs81srHPuSU+6815FPz/3AcOBpmbWzTn3ohfNBYCKfn6GUjBtGgnM8H9bAaXcMXLO3QtgZrcDe51z+f5urD4GvpVTc865I8Ad/m4mAFU0PvuAu/zdTACqaHyepeCgIdhVND7zgHn+bSVglTtGRX9xbrL/WimpPk5tZAAdi213AHZ61Esg0vhUTuNTOY3PqQXsGNXHwF8GdDezeDOLAG4G3ve4p0Ci8amcxqdyGp9TC9gxqtOBb2ZTgMVADzPLMLPRzrlc4F7gY2AdMM05t9bLPr2i8amcxqdyGp9Tq2tjpMXTRESCRJ0+whcRkapT4IuIBAkFvohIkFDgi4gECQW+iEiQUOCLiAQJBb6ISJBQ4IuIBAkFvohIkPh/DOg6wnH96ZMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cl_output_map = hp.anafast(output_map[0], lmax=16384)\n",
    "plt.loglog(cl_output_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "del map_log_pol_tens"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Galactic plane fix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "galplane_fix = hp.read_map(datadir / \"gnilc_dust_galplane.fits.gz\", (0, 1, 2, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_map *= hp.ud_grade(galplane_fix[3], output_nside)\n",
    "output_map += hp.ud_grade(galplane_fix[:3] * (1 - galplane_fix[3]), output_nside)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.72037508e+01,  3.66491306e+01,  3.67077472e+01, ...,\n",
       "         4.10425484e+01,  4.07681616e+01,  4.11603859e+01],\n",
       "       [ 7.03200141e-01, -6.66901679e-01,  6.07019814e-01, ...,\n",
       "         7.62817241e-01, -6.98765800e-01,  7.95575034e-01],\n",
       "       [-8.09556476e-02,  1.21250909e-01,  8.08351243e-03, ...,\n",
       "        -8.87983650e-01,  8.63559839e-01, -9.42672655e-01]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "output_map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Color correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Planck 353 GHz color correction https://github.com/galsci/pysm/issues/99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_map *= 0.911"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_map_sm = hp.smoothing(output_map, fwhm = np.radians(0.9/60))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7faae1fdf310>]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlcElEQVR4nO3deXxU1f3/8dcnu+xCWBSCCRg2sYgEZBFBBUGx0tp+FfRnXVDA1lZr1Qb3pQptv639Wm2FVoq1FYq41EoUFRWssgWUfYssJoAkCIQ1hJDz+yNxzE7CZObOZN7PxyOPx9wzd+Z+5hjf3Jw591xzziEiIg1flNcFiIhIcCjwRUQihAJfRCRCKPBFRCKEAl9EJEIo8EVEIkSM1wXUJDEx0SUnJ3tdhohIWFm+fPke51zriu0hHfjJyclkZmZ6XYaISFgxs+1VtWtIR0QkQijwRUQihAJfRCRCKPBFRCKEAl9EJEIo8EVEIoQCX0QkQijwRUQihAJfRCRCKPBFRCKEAl9EJEIo8EVEIoQCX0QkQijwRUQihAJfRCRCKPBFRCKEAl9EJEIELfDNbKiZfWxmz5vZ0GAdV0RESvgV+GY23cxyzWxNhfaRZrbRzLLMLL202QGHgAQgx5/jesE5x8rs/Rw6VuR1KSIip8Tfe9rOAJ4F/v5Ng5lFA88BwykJ9mVm9ibwsXNugZm1BX4PXO/nsYPi60PHeHVFDjOXZrN1z2HOS2rBzNv6c1pctNeliYjUiV9n+M65hcDeCs39gCzn3BbnXCEwCxjtnCsufX4fEO/PcQOtuNjxSdYefvLyCvpPns9TGRto1TiOn1zcmZU5+/nZrM84Uey8LlNEpE78PcOvSnsgu8x2DnCBmV0NjABaUPJXQZXMbDwwHqBjx44BKK96eQePMWd5DrOWfcn2r4/Q/LRYbuifzNh+SaS2bQpAm6YJPPLmWh77z1oeu+oczCyoNYqInKpABH5VCeicc68Br53sxc65acA0gLS0tICfRhcXOz75Yg8zl37Ju2t3U1Ts6JfSkp8P68LInu1IiC0/dHPjwGR27D/KtIVbaN/iNCYM6RzoEkVE6kUgAj8HSCqz3QHYGYDj+CX3QAGvlJ7NZ+89yumNYrlpYDJj+nXk7DZNanxt+shu7Nx/lMlvb6Bd8wRGn9c+SFWLiJy6QAT+MiDVzFKAHcAY4LoAHKfOThQ7Pt6cx8ylX/L++lxOFDsGdGrFPZd1ZcQ5lc/mqxMVZfzuml7kHjzGva+sok3TBAZ0bnVKNR0sOM6bK3fy/d7taRQXiP8cIiIl/EoYM5sJDAUSzSwHeMQ594KZ3QHMA6KB6c65tX5X6oev8gt4JTObWcuy2bH/KC0bx3HrhSlc2zeJTq1rPpuvTnxMNH+5IY0fPP8p41/K5NXbB9KldJy/tvKPHOdHf1vKyuz9bMk7zENX9jilWkREasOcC93ZJmlpaS4zM/OUXnui2LFgUy4vL8nmw40lZ/ODzm7F2H4dGd6jLfEx9TOtMmffEb7/p0+JjTJe/8kg2jZLqNXrvj50jBteWEpW7iHOS2rB8i/3kfGzwXRtV7d/NEREKjKz5c65tErtoRz4bTr1cNc++c86v845x4rt+9iZX0Bikzh+2CeJMX2TSE5sHIAqYc2OfK6duoiOrRoze0J/mibE1rj/7gMFXP/XJeTsO8JffpTGOWc255LffUSXtk351/j+mvkjIn4Jy8BvltTVpd059ZRem9SyEdf2TWJY97bExQR+BYkFm/K4ZcYyBnZuxfSb+hIbXfUxc/Yd4fq/LmHPwWNMv6kvF3QqGft/ecmX3P/6av5w7Xl8r3fNXwJn5R5kxZf7uSYtqcb9RCQyhWXg+zOk44XZmdncN2cVP+zTgd/+8DuVztS37TnM9X9dwsGC47x4Sz96dzzd99yJYsfVf/qEnfkFfPCLIdX+lbD/SCGjnvkvO/YfZc7EAaQltwzoZxKR8FNd4Gu1zHp0TVoSd16aypzlOTz9/uZyz2XlHuSaqYs4UljEy7f1Lxf2ANFRxuOje7Ln0DH+UOG133DOce+cVeQeLKD5abH83/yq9xMRqYoCv57dNSyV/+nTgWfmb+Zfy74EYN3OA1w7dTEO+NeEAfRs37zK1/ZKasHYfh2Z8ek2Nnx1oNLzMz7dxnvrdpN+eXduH9qZjzfvYcWX+6p8r3fW7OLWF5dRdKK4yudFJPIo8OuZmfHU1edyUZfW3P/6GqYu+IIx0xYRHxPF7AkDTjp1897LutI0IYaH31hL2eG2VTn7eSpjPcO6t+WWQcnc0P8sTm8UyzNVnOXnHz3OA6+v4f31ucxbu7veP6OIhCcFfgDERkfxp+vPp2vbpkx+ewMtGsXxrwkDSKnFLKHTG8fxy5HdWLptL//+vOQC5QMFx7nj5c9IbBLv+26gcXwMtw7uxEcb81iZvb/cezz93ib2HSkksUkcMz7dWukYew8X8mnWnnr5rCISPhT4AdIkPoYZN/dlwkWdmD1hAEktG9X6tdemJdErqQVPZqznQMFxJr22mh37j/LM2N6c3jjOt9+PBpxF89Ni+eMH357lb/jqAC8t3s51F3Rk4pDOLNu2jzU78su9/y9fXcV1f13ClrxD/n9QEQkbCvwAatMsgUlXdKdd89pdjPWNqCjjidHnsOfQMa6dupi5q3Zx9/Au9K0wI6dpQizjLkzh/fW5rNmRj3OOR/69lmYJMdxzWVf+Jy2JRnHRzPh0m+81mdv28t66kmGevy/aXu798o8e52DB8VP7sCIS8hT4Ieo7HVpwXb+OrN91gMGpidxezaqcNw1KpmlCDM/M38x/Vu1iyda93DOiKy0axdH8tFh+cH4H3vx8J3sOHcM5x5S3N9C6aTyX9WjLnOU5vjt4HT9RzA/+/CnjXgyfabAiUjcK/BB238hu3DUslT9cex5RUVVffdssIZZbBqXw7rrdPPLvNfRs34wxfb+9j8CNA8+i8EQxs0oXjMvcvo+7hqUycWhnDh0r4rUVJXebfGnRdrJyD7F0615W5ewHSqaBPv3eJhZv+Trgn1VEAk+BH8KanxbLXcO60KpJzTcIu2VQCk3jY9h35DiPXdWT6DL/OJzdpimDUxN5afF2fvPOBjolNuaatCR6J7WgV4fmvPjpNvYeLuQP72+iX0pLmsTH8LdPtgHw0cY8/m/+Zp6cu973fkcKizhaeCIgn1dEAkuB3wA0bxTLr3/4HR75bg/6nHV6pedvGZTC7gPH2Jx7iHtHdCU2Ogoz48aByXyRd5ibZyzj0LEinhjdkx/26cBbq3ayK/8oU97eQHSUsXpHPiuz91Nc7BgzbTE3z1jqwacUEX8p8BuIK849g5sHpVT53JAurencujG9O7ZgZM92vvZR3zmDxCZxrMzez3UXdKRru6bcNDCZomLHLTMy2bj7IE99vyeN46J5afF2/r1yB6ty8lm8ZS+bdx/kSGERN05fyqIvNOQjEg50x40IEBVlzJk4kJhoK7e+T3xMNDcPSmH6f7fy82FdAEhObMyl3drw/vpceiW14Jq0JFbl5DNneQ6fZu3h7DZN2P71YWYty6ZN03gWbMojJsoY0LkVi774mijDtyCciIQWneFHiNMbx1W5INuPh3Zm0aRLy31PMGFIZ+JjonhwVHfMjP/X/yyOFRWzM7+AR77bg+E92vLaihyeX/AF0VHGgk157Nh/lNv/uZxfvLIS5xyHjhVxpLAomB9RRE5CZ/gRzsyIiyk/A6hvckvWPjaCmNIlnruf0YzBqYnEx0QxOLU1xQ4yVn8FwOSrz2XSa6uZ8FIm+48cZ/+R46zekc9Db6yhSUIM/7y1f9A/k4hUTYEvVYqpsJ7/izf38z0efHYiKYmNObtNE8b0TeJvn2xlzY4D9GzfjA27DvLEW+tYmVNyde+aHfk88Ppqfnzx2Yw4px0i4h0N6UitREWZ71qAqCjj33cM4o9je2Nmvhu23D28CwM6t2LZtn00S4jBDMa9uIyVOfm88N+tvLzkS67+0yecKA7dezCINGRBO8M3s+7AnUAiMN859+dgHVvqX7My3wfcMiiF1DZNubhrG3YfOMbHm/dw48BkPs/ez8eb99A4LpqlW/ey8auD5B89zqvLc/jjh5v58/V9ql0qWkTqn19n+GY23cxyzWxNhfaRZrbRzLLMLB3AObfeOTcRuAaodCcWCV8JsdEM79EWM+OqXmcy4aJOjLswhR8NSKZpfAzP39AHKFmrJ8rgoX+vIXvvUf65ZDtjpi3iww25Hn8Ckcjg1y0Ozewi4BDwd+dcz9K2aGATMBzIAZYBY51z68zsKiAdeNY59/LJ3j/cbnEolTnnMDNueGEJJ4odjeJieH99+TX6z+/YgqPHi5lwUSe+2+vMclcKi0jdBeyetmaWDLxVJvAHAI8650aUbk8CcM5NLvOauc65UdW833hgPEDHjh37bN++vardJMwUFhXjcLy3bjc/m/kZNw5M9i3h8I2UxMbkHTzGvSO6snXPYe4b2ZVGcZpXIFJX1QV+IP5vag9kl9nOAS4ws6HA1UA8kFHdi51z04BpUHKGH4D6xANxMSWjh6POPYN+KS2Jj4nmk6w9jOnbkcffWkd8TBRb9xwG4JE31wKwftcBurZrysNX9iA6yirdFF5E6iYQgV/V/5XOOfcR8FEAjidhxMxo07Tk/gDv/nwIAK2axHFmi9O4duoimibEkn+0ZE3+JVv3smTrXv6+aDsPXNGdkT3b1elGMiJSXiCmZeYASWW2OwA7A3AcaSBGn9eevskteeung5kzcQBnNE+gbbPyK4Q+mbGewb/5kDc+2+H7B0FE6iYQZ/jLgFQzSwF2AGOA6wJwHGlgepzZDIBFky4l7+Axlmz9mllLs/lvmfvv3vWvzwEYf1En0kd2q/Y+ASJSmb+zdGYCQymZW78beMQ594KZXQH8AYgGpjvnnjyV99csHQEoOH6CfyzezvpdB3m19IYt37jz0lRuHpRMi0Zx1bxaJPIEbJZOICnwpSLnHF/kHeahN9awqPROXK2bxvPId3twabe2nBYX7XGFIt5T4EuDs37XAab/dyv/WbWTguPFANw+tDO/GN6l0lpAIpFEgS8N1oGC4/xyzireXvOVr+1nl6by46GdSYjVGb9EHgW+NHi7DxSQsXoXb3y2g5U5+TSJj2FIl9bcNSyV1LZNvS5PJGgU+BJRPt6cx+srdvDe+t0UFhUz+rwzuXVwJ7oo+CUCKPAlIuUeLOChN9Ywb23J+j3xMVHMuLkfAzrrNozScFUX+PpmSxq0Nk0TmHpDGq9MHMDpjWI5VlTM2L8spvtD7/DWqp2E8gmPSH3TGb5ElPnrdzPuxW9/p/olt+TmQckM69GWWM3skQZCQzoipZxzrNt1gH8s3s7rn+3wTen8/TW9+N557XX1roQ9Bb5IFfYdLuSeV1Yyv/QmLOd3bMHdw7tyYWqix5WJnDoFvkgN1u7M58E31rB+1wHfGf8LN6Zxafe2HlcmUncKfJFayD1QwPV/XcLm3EO+tiX3X0rbZgkeViVSN5qlI1ILbZol8N7dQ5haeh9egAuems87Za7iFQlXCnyRKow4px1bJ1/B46PPoVlCDBP/sZw7Z33GwQKtxS/hS4EvUg0z40cDklly/zBuG5zCW6t2ce6j7/Lch1lelyZyShT4IidxWlw0D4zqwazx/QH47byN3PpiJrvyj3pcmUjdKPBFaqlvcks2/mokE4d0ZuHmPC753wU88dY6Co6f8Lo0kVpR4IvUQXxMNOmXd2P+3UO4MDWRF/67lW4PvcMHG3Z7XZrISSnwRU5BUstGTLuhDwNLF2G7ZUYm976y0uOqRGoWiJuYV8nMOgEPAM2dcz8M1nFFAsXMePm2/uzcf5SfzfyMV5bnsP6rA1zX7yyuu6Cj1+WJVOLXGb6ZTTezXDNbU6F9pJltNLMsM0sHcM5tcc6N8+d4IqHozBanMXN8f269MIU1Ow5w/+urSU6fy/avD3tdmkg5/g7pzABGlm0ws2jgOeByoAcw1sx6+HkckZAWGx3Fg1f24KN7hvrahvz2IxZsyvOuKJEK/Ap859xCYG+F5n5AVukZfSEwCxhd2/c0s/FmlmlmmXl5+p9FwktyYmO2TRnFiHNK1uC5cfpSLnt6ASeKQ3cJE4kcgfjStj2QXWY7B2hvZq3M7Hmgt5lNqu7Fzrlpzrk051xa69atA1CeSOBNvSGN7/duD8Cm3YfofH8GS7dWPDcSCa5ABH5Vi4k759zXzrmJzrnOzrnJATiuSEh5+trz2DZlFBOGdALgmqmLmLn0S4+rkkgWiMDPAZLKbHcAdgbgOCJhYdLl3blrWGrJ49dW87/zNurWiuKJQAT+MiDVzFLMLA4YA7wZgOOIhI27hnXhlYkDiIuJ4tkPs0iZlMHRQl2hK8Hl77TMmcAioKuZ5ZjZOOdcEXAHMA9YD8x2zq31v1SR8NY3uSWfpl/i2+7+8Dt8vFkTEyR4dAMUEQ8892EWv523EYCfXXI2d1/W1eOKpCHRDVBEQshPLj6bP1x7HgDPfJDFaytyvC1IIoICX8Qj3+vdnoX3XgzA3bNXkpw+l/yjusGKBI4CX8RDHVs14vOHh/u2ez32Ljv2a519CQwFvojHWjSKY8MT365QMnbaYlbl7PeuIGmwFPgiISAhNpptU0bxt5v68uXeI1z17Ce8tHi712VJA6PAFwkhF3drw8NXlqw1+NAba/jj/M0eVyQNiQJfJMTccmEKiyddCsDv3tvE2fdnkHfwmMdVSUOgwBcJQe2aJ/DGTwYBUFTs6Pvk+xQWFXtclYQ7Bb5IiDovqQVbnrqCHmc0A6DLg29r2qb4RYEvEsKiooyMOwf71tfv9di7fJVf4HFVEq4U+CJhYOoNaXRt2xSA/pPns2n3QY8rknCkwBcJE/N+fpHv8WVPL+Tz7P3eFSNhSYEvEka2TRlFn7NOB+B7z33CQ2+s8bgiCScKfJEw8+rtA/nB+R0AeGnxdg3vSK0p8EXC0O+u6UWT+BigZHgn96C+yJWTU+CLhKnVj17me9zvyfkeViLhQoEvEqbMjK2Tr/BtJ6fP5WCB5ulL9RT4ImHMzNj4q29X2jz30Xc1pi/VClrgm1knM3vBzOYE65gikSA+Jrrcmf5lTy+kuDh0b10q3qlV4JvZdDPLNbM1FdpHmtlGM8sys/Sa3sM5t8U5N86fYkWkamZWbk39TvdncPhYkYcVSSiq7Rn+DGBk2QYziwaeAy4HegBjzayHmZ1rZm9V+GlTr1WLSCUJsSVn+n2TS+bpn/PIPIW+lFOrwHfOLQT2VmjuB2SVnrkXArOA0c651c65Kyv85NZz3SJSBTNj9oQBvu1zHpnnYTUSavwZw28PZJfZziltq5KZtTKz54HeZjaphv3Gm1mmmWXm5eX5UZ5IZKo4e6fLg297WI2EEn8C36poq/abIufc1865ic65zs65yTXsN805l+acS2vdurUf5YlELjPjrZ9eCEBhUTGP/2edxxVJKPAn8HOApDLbHYCd/pUjIvWlZ/vm/PPWCwCY/slWRj3zsccVidf8CfxlQKqZpZhZHDAGeLN+yhKR+jDo7ERe//FAANbuPEBy+lyPKxIv1XZa5kxgEdDVzHLMbJxzrgi4A5gHrAdmO+fWBq5UETkVvTuezm2DU3zbGat3eViNeMmcC90LNNLS0lxmZqbXZYg0CBNfWs47a78CYOvkKzCr6ms4aQjMbLlzLq1iu5ZWEIkQz9/Qh+/3LplIlzIpw+NqxAsKfJEI8qvv9fQ91nh+5FHgi0SQxvExvHDjt3/p/++8jR5WI8GmwBeJMJd2b8u9I7oC8OyHWYTy93hSvxT4IhHoJxef7Xt8zyurPKxEgkmBLxKh1j0+AoBXV+RoPD9CKPBFIlSjuBh+cnFn33b23iMeViPBoMAXiWD3juhGp9aNARj8mw89rkYCTYEvEuE++MVQ32MN7TRsCnwR4bnrzvc9Vug3XAp8EWHUd84ot52zT+P5DZECX0QA2DZllO/xhb/+kE27D3pYjQSCAl9EfLKevNz3+LKnF7L3cKGH1Uh9U+CLiE9MdBSrH73Mt33+E+9xtPCEhxVJfVLgi0g5TRNifbdHBOj+8DvkHz3uYUVSXxT4IlJJz/bN+eieob7tXo+9y39W6g6m4U6BLyJVSk5szIYnRvq2fzrzMwZN+UCLrYUxBb6IVCshNprNZb7I3bH/KCmTMig4rnH9cKTAF5EaxUZH+RZa+0a3h94hOX2uxvbDTNAC38y6m9nzZjbHzG4P1nFFxH+N4mLYNmUUM2/rX66912PvsuLLfRrmCRO1Cnwzm25muWa2pkL7SDPbaGZZZpZe03s459Y75yYC1wCVbq4rIqFvQOdWvHPX4HJtV//pU1ImZbB1z2GPqpLastr8y2xmFwGHgL8753qWtkUDm4DhQA6wDBgLRAOTK7zFLc65XDO7CkgHnnXOvXyy46alpbnMzMw6fBwRCZbjJ4pJfeDtSu1fPHUF0VHmQUXyDTNb7pyrdGJdq8AvfYNk4K0ygT8AeNQ5N6J0exKAc65i2Ff1XnOdc6OqeW48MB6gY8eOfbZv316r+kQk+AqOn+DOWZ8xb+3ucu09zmhGxp2Dq3mVBFp1ge/PGH57ILvMdk5pW3UFDDWzZ8xsKpBR3X7OuWnOuTTnXFrr1q39KE9EAi0hNpqpN6SVm74JsG7XAZLT5zJ7WXY1rxQvxPjx2qr+Zqv2zwXn3EfAR34cT0RCVEJsNNumjOJY0Qm6PviOr/2+V1dx36ur6JfSktkTBnhYoYB/Z/g5QFKZ7Q6ALsUTiWDxMSXBP/nqc8u1L926l+T0uRQXazaPl/wJ/GVAqpmlmFkcMAZ4s37KEpFwNrZfR7ZNGcWEizqVa+90f4ZusOKh2k7LnAksArqaWY6ZjXPOFQF3APOA9cBs59zawJUqIuFm0hXd2TZlFFf1OrNce3L6XB54fTV7Dh3zqLLIVOtZOl7QtEyRhuWypxewafehcm0TLupE+uXdMNNUzvoSiFk6IiJ18u7Ph1RapmHqwi2kTMrQbRWDQIEvIkH1zTINH993cbn2C3/9IZPfXu9RVZFBgS8inkhq2ajki90h336xO3XBFpLT53L1nz7R+jwBoDF8EfGcc46USZWvx0w763Tm3D7Qg4rCm8bwRSRkmRnbpoxi/ePlr9jN3L6P5PS5Gt+vJzrDF5GQc6LY0fn+ymf8WpitdnSGLyJhIzqq5Ix/6g19yrV3vj+DG6cv9aiq8KfAF5GQNeKcdmybMorhPdr62hZsyiM5fa7W3z8FGtIRkbCw73AhvZ94r1L7tilVrrQe0TSkIyJh7fTGcVXO309On0vvx9/1qKrwosAXkbDyzfz9Tq0b+9r2HTlOcvpcFm7K87Cy0KchHREJWwcLjnPuo5XP7rc8dQVRETybR0M6ItLgNE2IZduUUcz/xZBy7Z3uz+Ca5xd5VFXoUuCLSNjr3LoJ26aM4vKe7XxtS7eV3HRlgYZ5fDSkIyINSnGxo1MVF21tnXxFxCzBrCEdEYkIUaUXbb1/90Xl2lMmZXDPKys9qio0KPBFpEE6u03TSnP05yzPITl9LseKTnhUlbcU+CLSoG2bMootT11Rrq3rg+/wi9mRd7YftMA3s6Fm9rGZPW9mQ4N1XBGRb4Z5Bqcm+tpeXVFytn/8RLGHlQVXbW9iPt3Mcs1sTYX2kWa20cyyzCz9JG/jgENAApBzauWKiJy6l8ZdUGmYJ/WBt/nLwi0eVRRctZqlY2YXURLWf3fO9SxtiwY2AcMpCfBlwFggGphc4S1uAfY454rNrC3we+fc9Sc7rmbpiEigfLBhN7fMKJ8vDeWCLb9m6TjnFgJ7KzT3A7Kcc1ucc4XALGC0c261c+7KCj+5zrlv/m7aB8TXUOh4M8s0s8y8PM2fFZHAuKRb20pj+53uz+CfS7Z7VFHg+TOG3x7ILrOdU9pWJTO72symAi8Bz1a3n3NumnMuzTmX1rp1az/KExGp2Tdj+3++/nxf2wOvryE5fa6HVQWOP4Ff1d891Y4POedec85NcM5d65z7yI/jiojUq8vPPYMvKpztJ6fP5etDx+r9WCuz9/PBht31/r614U/g5wBJZbY7ADv9K0dExBvf3GUrsUmcr63Pr97npzM/q9fjjH7uE26ZkenJDVz8CfxlQKqZpZhZHDAGeLN+yhIR8Ubmg8P5/OHhvu3/rNxJcvpc6nsZmjnLs0++Uz2r7bTMmcAioKuZ5ZjZOOdcEXAHMA9YD8x2zq0NXKkiIsHRolFcpembKZMyyD9yvN6O8fWhwnp7r9qq7Sydsc65M5xzsc65Ds65F0rbM5xzXZxznZ1zTwa2VBGR4No2ZRQPXdnDt93r8Xd5JbN+zsyPFQX/gi8trSAiUoNxF6aw+tHLfNv3zll1yrN4Co5/u4bP0cLgr+ejwBcROYmmCbFsnVx5Fk9d7S8zJHS4sMjvuupKgS8iUgtmVmlcv66hvzP/qO/xwQIFvohISNs2ZRSTLu/m205On8vBgtp9mVt2SKe2r6lPCnwRkTqaMKQzmQ8O822f++i77Cpz9l6dA0e/Dfkv8sJrHr6ISMRKbBLPusdH+LYHTP6ALXmHanzNgaPBH8YpS4EvInKKGsXFkPXk5b7tS363gIzVu6rdf9+R4M+9L0uBLyLih5joqHIzeH78zxW88dmOKvfde1iBLyIS1irO4LnrX5+z4st9lfarGPiFQb74SoEvIlJPyob+1X/6tNKY/hcVtoN9xq/AFxGpR2VD/5LfLSB77xHf9oov95fbd08All+uiQJfRKSelR3TH/ybDzlSzVW1RcX1uwLnySjwRUTqmZmVC/0eD8+rcnnlQ0G+2laBLyISAGbG5jJTNlMmZVTa50CQr7ZV4IuIBEhsdBQL7h1a7fNlr7wNBgW+iEgAndWqMU9+v6dvOz7m29itaupmIMUE9WgiIhHo+gvOYsQ57ThYUMSqnP3c+8oqCk8Us3F3zUsx1Ded4YuIBEFik3hSEhsz+rz2bCod2/86yNMyg3aGb2aDgetLj9nDOTcwWMcWEQlFOftOvsJmfartTcynm1muma2p0D7SzDaaWZaZpdf0Hs65j51zE4G3gBdPvWQRkfDXK6lF0I9Z2yGdGcDIsg1mFg08B1wO9ADGmlkPMzvXzN6q8NOmzEuvA2bWQ+0iImFrZfZ+ALJyDwbtmLUKfOfcQmBvheZ+QJZzbotzrhCYBYx2zq12zl1Z4ScXwMw6AvnOuQP1+SFERMLNj4d2BmDY7xcG7Zj+fGnbHsgus51T2laTccDfatrBzMabWaaZZebl5flRnohI6Lpv5Le3SRwweX5QjulP4FsVbTUuDOGce8Q59+lJ9pnmnEtzzqW1bt3aj/JERELbP8ZdAMCu/AIOHwv8Mgv+BH4OkFRmuwOw079yREQix4Wpib7H5zwyL+DH8yfwlwGpZpZiZnHAGODN+ilLRCQylF1k7RezVwb0WLWdljkTWAR0NbMcMxvnnCsC7gDmAeuB2c65tYErVUSk4TEz7hvZFYBXV+QE9lhVLdkZKtLS0lxmZqbXZYiIBFxy+lzf47I3UTkVZrbcOZdWsV1LK4iIhIDFky71PS44fiIgx1Dgi4iEgHbNE3yPuz30TkCOocAXEQkRZb/ADcRZvgJfRCREmH17edPanfW/IIECX0QkBB0t1Bm+iEhEOFxY/1feKvBFRELQEQW+iEhkOHRMQzoiIg3anZemAtCn4+n1/t660lZEpIHRlbYiIhFOgS8iEiEU+CIiEUKBLyISIRT4IiIRQoEvIhIhFPgiIhFCgS8iEiFC+sIrM8sDtpdpag7k13I7EdgTgLIqHrO+9q9pv6qeq01bOPRPbV+j/jn1fap7rq59VPG5QPRROPdPxW0vf4fOcs61rtTqnAubH2BabbeBzGDUUF/717RfVc/Vpi0c+qe2r1H/1G//nEofVfFcvfdROPdPKP8OffMTbkM6/6njdjBqqK/9a9qvqudq0xYO/VPb16h/Tn2f6p6rax+pf8L3dwgI8SEdf5hZpqtiLQkpof6pmfrn5NRHNQvF/gm3M/y6mOZ1ASFO/VMz9c/JqY9qFnL902DP8EVEpLyGfIYvIiJlKPBFRCKEAl9EJEJETOCbWWMze9HM/mJm13tdT6gxs05m9oKZzfG6llBkZt8r/d35t5ld5nU9ocbMupvZ82Y2x8xu97qeUFSaQcvN7EqvagjrwDez6WaWa2ZrKrSPNLONZpZlZumlzVcDc5xztwFXBb1YD9Slf5xzW5xz47yp1Bt17J83Sn93bgKu9aDcoKtj/6x3zk0ErgFCaipioNQxfwB+CcwObpXlhXXgAzOAkWUbzCwaeA64HOgBjDWzHkAHILt0t/q/HXxomkHt+ycSzaDu/fNg6fORYAZ16B8zuwr4LzA/uGV6Zga17B8zGwasA3YHu8iywjrwnXMLgb0VmvsBWaVnrIXALGA0kENJ6EOYf+7aqmP/RJy69I+V+DXwtnNuRbBr9UJdf3+cc2865wYCETFkWsf+uRjoD1wH3GZmnmRQjBcHDbD2fHsmDyVBfwHwDPCsmY0iOJc/h6oq+8fMWgFPAr3NbJJzbrIn1Xmvut+fnwLDgOZmdrZz7nkvigsB1f3+DKVk2DQeyAh+WSGjyv5xzt0BYGY3AXucc8Ue1NYgA9+qaHPOucPAzcEuJgRV1z9fAxODXUwIqq5/nqHkpCHSVdc/HwEfBbeUkFRl//geODcjeKVU1hCHNnKApDLbHYCdHtUSitQ/NVP/1Ez9U7OQ7p+GGPjLgFQzSzGzOGAM8KbHNYUS9U/N1D81U//ULKT7J6wD38xmAouArmaWY2bjnHNFwB3APGA9MNs5t9bLOr2i/qmZ+qdm6p+ahWP/aPE0EZEIEdZn+CIiUnsKfBGRCKHAFxGJEAp8EZEIocAXEYkQCnwRkQihwBcRiRAKfBGRCKHAFxGJEP8fZPW2XinZ048AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cl_output_map_sm = hp.anafast(output_map_sm[0], lmax=16384)\n",
    "plt.loglog(cl_output_map_sm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl_output_map_postfix = hp.anafast(output_map[0], lmax=16384)\n",
    "plt.loglog(cl_output_map_postfix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(galplane_fix[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_map[0][hp.ud_grade(galplane_fix[3], 8192)<1] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(output_map[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl_output_map_postfix_masked = hp.anafast(output_map[0], lmax=16384)\n",
    "plt.loglog(cl_output_map_postfix_masked)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.write_map(\n",
    "    output_dir / f\"gnilc_dust_template_nside{output_nside}.fits\",\n",
    "    output_map,\n",
    "    dtype=np.float32,\n",
    "    overwrite=True,\n",
    "    column_units = \"uK_RJ\",\n",
    "    extra_header = [(\"lmax\", output_lmax), (\"ref_freq\", \"353 GHz\")]\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "condanamaster2",
   "language": "python",
   "name": "condanamaster2"
  },
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}