FIRAS_exercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "11280ebb",
   "metadata": {},
   "source": [
    "<div align=\"center\"><h1>\n",
    " AST9240/AST5240 Project __\n",
    "    </br>\n",
    "   FIRAS: Sampling for spectral distortions\n",
    "</h1></div>\n",
    "<hr style=\"border:2px solid gray\"> </hr>\n",
    "\n",
    "#### Header and imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f955f114",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import emcee as mc\n",
    "from corner import corner\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy import constants as const\n",
    "from astropy.table import QTable\n",
    "from astropy.modeling.models import BlackBody\n",
    "\n",
    "from IPython.display import display, Math\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib_inline.backend_inline import set_matplotlib_formats\n",
    "from astropy.visualization import quantity_support\n",
    "\n",
    "quantity_support()\n",
    "set_matplotlib_formats('svg')\n",
    "from matplotlib import cm\n",
    "from cycler import cycler\n",
    "plt.rc('legend', frameon = False)\n",
    "plt.rc('figure', figsize = (7, 7/1.25))\n",
    "plt.rc('font', size = 14)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10936c4e",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "We're using data from Table 4 of [Fixsen et. al. (1996)](https://articles.adsabs.harvard.edu/pdf/1996ApJ...473..576F), which consists of frequencies in cm$^{-1}$, as well as residuals, uncertainties, and a model of the Milky Way, all in kJy/sr. There is also an artificial monopole spectrum included, which is produced by adding a 2.725 K blackbody spectrum to the given residuals. This seems to be the current accepted value for T$_0$, although it is not the one actually reported in Fixsen et. al. (1996). \n",
    "\n",
    "It might be useful to make an alternate array of frequencies in s$^{-1}$, depending on how you plan on formatting your models, but is not technically necessary."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a443d4c7-9c59-419d-8c67-d870468cf70f",
   "metadata": {},
   "source": [
    "First, let's read in the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9777ef7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get the data\n",
    "names = [\"freq\", \"monopole\", \"res\", \"sigma\", \"gal\"]\n",
    "tmp = pd.read_csv(\"https://lambda.gsfc.nasa.gov/data/cobe/firas/monopole_spec/firas_monopole_spec_v1.txt\", delim_whitespace = True, header = 0, skiprows = 17, names = names)\n",
    "unit_names = [\"cm-1\", \"MJy/sr\", \"kJy/sr\", \"kJy/sr\", \"kJy/sr\"]\n",
    "units = {names[i]:unit_names[i] for i in range(len(names))}\n",
    "firas = QTable.from_pandas(tmp, units = units)\n",
    "\n",
    "#make a version of the frequency in s-1\n",
    "firas_freq = firas[\"freq\"].to('GHz', equivalencies = u.spectral())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f7873f5-5158-44d0-8c7b-1221eed6c58c",
   "metadata": {},
   "source": [
    "and plot the monopole data, with the quoted $T_0$ blackbody to compare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9a239b87",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       "  \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<svg xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"456.271875pt\" height=\"375.680125pt\" viewBox=\"0 0 456.271875 375.680125\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">\n",
       " <metadata>\n",
       "  <rdf:RDF xmlns:dc=\"http://purl.org/dc/elements/1.1/\" xmlns:cc=\"http://creativecommons.org/ns#\" xmlns:rdf=\"http://www.w3.org/1999/02/22-rdf-syntax-ns#\">\n",
       "   <cc:Work>\n",
       "    <dc:type rdf:resource=\"http://purl.org/dc/dcmitype/StillImage\"/>\n",
       "    <dc:date>2022-08-30T13:50:59.049561</dc:date>\n",
       "    <dc:format>image/svg+xml</dc:format>\n",
       "    <dc:creator>\n",
       "     <cc:Agent>\n",
       "      <dc:title>Matplotlib v3.5.1, https://matplotlib.org/</dc:title>\n",
       "     </cc:Agent>\n",
       "    </dc:creator>\n",
       "   </cc:Work>\n",
       "  </rdf:RDF>\n",
       " </metadata>\n",
       " <defs>\n",
       "  <style type=\"text/css\">*{stroke-linejoin: round; stroke-linecap: butt}</style>\n",
       " </defs>\n",
       " <g id=\"figure_1\">\n",
       "  <g id=\"patch_1\">\n",
       "   <path d=\"M 0 375.680125 \n",
       "L 456.271875 375.680125 \n",
       "L 456.271875 0 \n",
       "L 0 0 \n",
       "L 0 375.680125 \n",
       "z\n",
       "\" style=\"fill: none\"/>\n",
       "  </g>\n",
       "  <g id=\"axes_1\">\n",
       "   <g id=\"patch_2\">\n",
       "    <path d=\"M 58.471875 330.381375 \n",
       "L 449.071875 330.381375 \n",
       "L 449.071875 25.965375 \n",
       "L 58.471875 25.965375 \n",
       "z\n",
       "\" style=\"fill: #ffffff\"/>\n",
       "   </g>\n",
       "   <g id=\"matplotlib.axis_1\">\n",
       "    <g id=\"xtick_1\">\n",
       "     <g id=\"line2d_1\">\n",
       "      <defs>\n",
       "       <path id=\"mb047b0e9d3\" d=\"M 0 0 \n",
       "L 0 3.5 \n",
       "\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </defs>\n",
       "      <g>\n",
       "       <use xlink:href=\"#mb047b0e9d3\" x=\"96.079485\" y=\"330.381375\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_1\">\n",
       "      <!-- 100 -->\n",
       "      <g transform=\"translate(82.718235 348.019187)scale(0.14 -0.14)\">\n",
       "       <defs>\n",
       "        <path id=\"DejaVuSans-31\" d=\"M 794 531 \n",
       "L 1825 531 \n",
       "L 1825 4091 \n",
       "L 703 3866 \n",
       "L 703 4441 \n",
       "L 1819 4666 \n",
       "L 2450 4666 \n",
       "L 2450 531 \n",
       "L 3481 531 \n",
       "L 3481 0 \n",
       "L 794 0 \n",
       "L 794 531 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "        <path id=\"DejaVuSans-30\" d=\"M 2034 4250 \n",
       "Q 1547 4250 1301 3770 \n",
       "Q 1056 3291 1056 2328 \n",
       "Q 1056 1369 1301 889 \n",
       "Q 1547 409 2034 409 \n",
       "Q 2525 409 2770 889 \n",
       "Q 3016 1369 3016 2328 \n",
       "Q 3016 3291 2770 3770 \n",
       "Q 2525 4250 2034 4250 \n",
       "z\n",
       "M 2034 4750 \n",
       "Q 2819 4750 3233 4129 \n",
       "Q 3647 3509 3647 2328 \n",
       "Q 3647 1150 3233 529 \n",
       "Q 2819 -91 2034 -91 \n",
       "Q 1250 -91 836 529 \n",
       "Q 422 1150 422 2328 \n",
       "Q 422 3509 836 4129 \n",
       "Q 1250 4750 2034 4750 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       </defs>\n",
       "       <use xlink:href=\"#DejaVuSans-31\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"xtick_2\">\n",
       "     <g id=\"line2d_2\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#mb047b0e9d3\" x=\"158.22302\" y=\"330.381375\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_2\">\n",
       "      <!-- 200 -->\n",
       "      <g transform=\"translate(144.86177 348.019187)scale(0.14 -0.14)\">\n",
       "       <defs>\n",
       "        <path id=\"DejaVuSans-32\" d=\"M 1228 531 \n",
       "L 3431 531 \n",
       "L 3431 0 \n",
       "L 469 0 \n",
       "L 469 531 \n",
       "Q 828 903 1448 1529 \n",
       "Q 2069 2156 2228 2338 \n",
       "Q 2531 2678 2651 2914 \n",
       "Q 2772 3150 2772 3378 \n",
       "Q 2772 3750 2511 3984 \n",
       "Q 2250 4219 1831 4219 \n",
       "Q 1534 4219 1204 4116 \n",
       "Q 875 4013 500 3803 \n",
       "L 500 4441 \n",
       "Q 881 4594 1212 4672 \n",
       "Q 1544 4750 1819 4750 \n",
       "Q 2544 4750 2975 4387 \n",
       "Q 3406 4025 3406 3419 \n",
       "Q 3406 3131 3298 2873 \n",
       "Q 3191 2616 2906 2266 \n",
       "Q 2828 2175 2409 1742 \n",
       "Q 1991 1309 1228 531 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       </defs>\n",
       "       <use xlink:href=\"#DejaVuSans-32\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"xtick_3\">\n",
       "     <g id=\"line2d_3\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#mb047b0e9d3\" x=\"220.366555\" y=\"330.381375\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_3\">\n",
       "      <!-- 300 -->\n",
       "      <g transform=\"translate(207.005305 348.019187)scale(0.14 -0.14)\">\n",
       "       <defs>\n",
       "        <path id=\"DejaVuSans-33\" d=\"M 2597 2516 \n",
       "Q 3050 2419 3304 2112 \n",
       "Q 3559 1806 3559 1356 \n",
       "Q 3559 666 3084 287 \n",
       "Q 2609 -91 1734 -91 \n",
       "Q 1441 -91 1130 -33 \n",
       "Q 819 25 488 141 \n",
       "L 488 750 \n",
       "Q 750 597 1062 519 \n",
       "Q 1375 441 1716 441 \n",
       "Q 2309 441 2620 675 \n",
       "Q 2931 909 2931 1356 \n",
       "Q 2931 1769 2642 2001 \n",
       "Q 2353 2234 1838 2234 \n",
       "L 1294 2234 \n",
       "L 1294 2753 \n",
       "L 1863 2753 \n",
       "Q 2328 2753 2575 2939 \n",
       "Q 2822 3125 2822 3475 \n",
       "Q 2822 3834 2567 4026 \n",
       "Q 2313 4219 1838 4219 \n",
       "Q 1578 4219 1281 4162 \n",
       "Q 984 4106 628 3988 \n",
       "L 628 4550 \n",
       "Q 988 4650 1302 4700 \n",
       "Q 1616 4750 1894 4750 \n",
       "Q 2613 4750 3031 4423 \n",
       "Q 3450 4097 3450 3541 \n",
       "Q 3450 3153 3228 2886 \n",
       "Q 3006 2619 2597 2516 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       </defs>\n",
       "       <use xlink:href=\"#DejaVuSans-33\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"xtick_4\">\n",
       "     <g id=\"line2d_4\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#mb047b0e9d3\" x=\"282.51009\" y=\"330.381375\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_4\">\n",
       "      <!-- 400 -->\n",
       "      <g transform=\"translate(269.14884 348.019187)scale(0.14 -0.14)\">\n",
       "       <defs>\n",
       "        <path id=\"DejaVuSans-34\" d=\"M 2419 4116 \n",
       "L 825 1625 \n",
       "L 2419 1625 \n",
       "L 2419 4116 \n",
       "z\n",
       "M 2253 4666 \n",
       "L 3047 4666 \n",
       "L 3047 1625 \n",
       "L 3713 1625 \n",
       "L 3713 1100 \n",
       "L 3047 1100 \n",
       "L 3047 0 \n",
       "L 2419 0 \n",
       "L 2419 1100 \n",
       "L 313 1100 \n",
       "L 313 1709 \n",
       "L 2253 4666 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       </defs>\n",
       "       <use xlink:href=\"#DejaVuSans-34\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"xtick_5\">\n",
       "     <g id=\"line2d_5\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#mb047b0e9d3\" x=\"344.653625\" y=\"330.381375\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_5\">\n",
       "      <!-- 500 -->\n",
       "      <g transform=\"translate(331.292375 348.019187)scale(0.14 -0.14)\">\n",
       "       <defs>\n",
       "        <path id=\"DejaVuSans-35\" d=\"M 691 4666 \n",
       "L 3169 4666 \n",
       "L 3169 4134 \n",
       "L 1269 4134 \n",
       "L 1269 2991 \n",
       "Q 1406 3038 1543 3061 \n",
       "Q 1681 3084 1819 3084 \n",
       "Q 2600 3084 3056 2656 \n",
       "Q 3513 2228 3513 1497 \n",
       "Q 3513 744 3044 326 \n",
       "Q 2575 -91 1722 -91 \n",
       "Q 1428 -91 1123 -41 \n",
       "Q 819 9 494 109 \n",
       "L 494 744 \n",
       "Q 775 591 1075 516 \n",
       "Q 1375 441 1709 441 \n",
       "Q 2250 441 2565 725 \n",
       "Q 2881 1009 2881 1497 \n",
       "Q 2881 1984 2565 2268 \n",
       "Q 2250 2553 1709 2553 \n",
       "Q 1456 2553 1204 2497 \n",
       "Q 953 2441 691 2322 \n",
       "L 691 4666 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       </defs>\n",
       "       <use xlink:href=\"#DejaVuSans-35\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"xtick_6\">\n",
       "     <g id=\"line2d_6\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#mb047b0e9d3\" x=\"406.79716\" y=\"330.381375\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_6\">\n",
       "      <!-- 600 -->\n",
       "      <g transform=\"translate(393.43591 348.019187)scale(0.14 -0.14)\">\n",
       "       <defs>\n",
       "        <path id=\"DejaVuSans-36\" d=\"M 2113 2584 \n",
       "Q 1688 2584 1439 2293 \n",
       "Q 1191 2003 1191 1497 \n",
       "Q 1191 994 1439 701 \n",
       "Q 1688 409 2113 409 \n",
       "Q 2538 409 2786 701 \n",
       "Q 3034 994 3034 1497 \n",
       "Q 3034 2003 2786 2293 \n",
       "Q 2538 2584 2113 2584 \n",
       "z\n",
       "M 3366 4563 \n",
       "L 3366 3988 \n",
       "Q 3128 4100 2886 4159 \n",
       "Q 2644 4219 2406 4219 \n",
       "Q 1781 4219 1451 3797 \n",
       "Q 1122 3375 1075 2522 \n",
       "Q 1259 2794 1537 2939 \n",
       "Q 1816 3084 2150 3084 \n",
       "Q 2853 3084 3261 2657 \n",
       "Q 3669 2231 3669 1497 \n",
       "Q 3669 778 3244 343 \n",
       "Q 2819 -91 2113 -91 \n",
       "Q 1303 -91 875 529 \n",
       "Q 447 1150 447 2328 \n",
       "Q 447 3434 972 4092 \n",
       "Q 1497 4750 2381 4750 \n",
       "Q 2619 4750 2861 4703 \n",
       "Q 3103 4656 3366 4563 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       </defs>\n",
       "       <use xlink:href=\"#DejaVuSans-36\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"text_7\">\n",
       "     <!-- Frequency (GHz) -->\n",
       "     <g transform=\"translate(195.561406 365.568562)scale(0.14 -0.14)\">\n",
       "      <defs>\n",
       "       <path id=\"DejaVuSans-46\" d=\"M 628 4666 \n",
       "L 3309 4666 \n",
       "L 3309 4134 \n",
       "L 1259 4134 \n",
       "L 1259 2759 \n",
       "L 3109 2759 \n",
       "L 3109 2228 \n",
       "L 1259 2228 \n",
       "L 1259 0 \n",
       "L 628 0 \n",
       "L 628 4666 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-72\" d=\"M 2631 2963 \n",
       "Q 2534 3019 2420 3045 \n",
       "Q 2306 3072 2169 3072 \n",
       "Q 1681 3072 1420 2755 \n",
       "Q 1159 2438 1159 1844 \n",
       "L 1159 0 \n",
       "L 581 0 \n",
       "L 581 3500 \n",
       "L 1159 3500 \n",
       "L 1159 2956 \n",
       "Q 1341 3275 1631 3429 \n",
       "Q 1922 3584 2338 3584 \n",
       "Q 2397 3584 2469 3576 \n",
       "Q 2541 3569 2628 3553 \n",
       "L 2631 2963 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-65\" d=\"M 3597 1894 \n",
       "L 3597 1613 \n",
       "L 953 1613 \n",
       "Q 991 1019 1311 708 \n",
       "Q 1631 397 2203 397 \n",
       "Q 2534 397 2845 478 \n",
       "Q 3156 559 3463 722 \n",
       "L 3463 178 \n",
       "Q 3153 47 2828 -22 \n",
       "Q 2503 -91 2169 -91 \n",
       "Q 1331 -91 842 396 \n",
       "Q 353 884 353 1716 \n",
       "Q 353 2575 817 3079 \n",
       "Q 1281 3584 2069 3584 \n",
       "Q 2775 3584 3186 3129 \n",
       "Q 3597 2675 3597 1894 \n",
       "z\n",
       "M 3022 2063 \n",
       "Q 3016 2534 2758 2815 \n",
       "Q 2500 3097 2075 3097 \n",
       "Q 1594 3097 1305 2825 \n",
       "Q 1016 2553 972 2059 \n",
       "L 3022 2063 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-71\" d=\"M 947 1747 \n",
       "Q 947 1113 1208 752 \n",
       "Q 1469 391 1925 391 \n",
       "Q 2381 391 2643 752 \n",
       "Q 2906 1113 2906 1747 \n",
       "Q 2906 2381 2643 2742 \n",
       "Q 2381 3103 1925 3103 \n",
       "Q 1469 3103 1208 2742 \n",
       "Q 947 2381 947 1747 \n",
       "z\n",
       "M 2906 525 \n",
       "Q 2725 213 2448 61 \n",
       "Q 2172 -91 1784 -91 \n",
       "Q 1150 -91 751 415 \n",
       "Q 353 922 353 1747 \n",
       "Q 353 2572 751 3078 \n",
       "Q 1150 3584 1784 3584 \n",
       "Q 2172 3584 2448 3432 \n",
       "Q 2725 3281 2906 2969 \n",
       "L 2906 3500 \n",
       "L 3481 3500 \n",
       "L 3481 -1331 \n",
       "L 2906 -1331 \n",
       "L 2906 525 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-75\" d=\"M 544 1381 \n",
       "L 544 3500 \n",
       "L 1119 3500 \n",
       "L 1119 1403 \n",
       "Q 1119 906 1312 657 \n",
       "Q 1506 409 1894 409 \n",
       "Q 2359 409 2629 706 \n",
       "Q 2900 1003 2900 1516 \n",
       "L 2900 3500 \n",
       "L 3475 3500 \n",
       "L 3475 0 \n",
       "L 2900 0 \n",
       "L 2900 538 \n",
       "Q 2691 219 2414 64 \n",
       "Q 2138 -91 1772 -91 \n",
       "Q 1169 -91 856 284 \n",
       "Q 544 659 544 1381 \n",
       "z\n",
       "M 1991 3584 \n",
       "L 1991 3584 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-6e\" d=\"M 3513 2113 \n",
       "L 3513 0 \n",
       "L 2938 0 \n",
       "L 2938 2094 \n",
       "Q 2938 2591 2744 2837 \n",
       "Q 2550 3084 2163 3084 \n",
       "Q 1697 3084 1428 2787 \n",
       "Q 1159 2491 1159 1978 \n",
       "L 1159 0 \n",
       "L 581 0 \n",
       "L 581 3500 \n",
       "L 1159 3500 \n",
       "L 1159 2956 \n",
       "Q 1366 3272 1645 3428 \n",
       "Q 1925 3584 2291 3584 \n",
       "Q 2894 3584 3203 3211 \n",
       "Q 3513 2838 3513 2113 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-63\" d=\"M 3122 3366 \n",
       "L 3122 2828 \n",
       "Q 2878 2963 2633 3030 \n",
       "Q 2388 3097 2138 3097 \n",
       "Q 1578 3097 1268 2742 \n",
       "Q 959 2388 959 1747 \n",
       "Q 959 1106 1268 751 \n",
       "Q 1578 397 2138 397 \n",
       "Q 2388 397 2633 464 \n",
       "Q 2878 531 3122 666 \n",
       "L 3122 134 \n",
       "Q 2881 22 2623 -34 \n",
       "Q 2366 -91 2075 -91 \n",
       "Q 1284 -91 818 406 \n",
       "Q 353 903 353 1747 \n",
       "Q 353 2603 823 3093 \n",
       "Q 1294 3584 2113 3584 \n",
       "Q 2378 3584 2631 3529 \n",
       "Q 2884 3475 3122 3366 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-79\" d=\"M 2059 -325 \n",
       "Q 1816 -950 1584 -1140 \n",
       "Q 1353 -1331 966 -1331 \n",
       "L 506 -1331 \n",
       "L 506 -850 \n",
       "L 844 -850 \n",
       "Q 1081 -850 1212 -737 \n",
       "Q 1344 -625 1503 -206 \n",
       "L 1606 56 \n",
       "L 191 3500 \n",
       "L 800 3500 \n",
       "L 1894 763 \n",
       "L 2988 3500 \n",
       "L 3597 3500 \n",
       "L 2059 -325 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-20\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-28\" d=\"M 1984 4856 \n",
       "Q 1566 4138 1362 3434 \n",
       "Q 1159 2731 1159 2009 \n",
       "Q 1159 1288 1364 580 \n",
       "Q 1569 -128 1984 -844 \n",
       "L 1484 -844 \n",
       "Q 1016 -109 783 600 \n",
       "Q 550 1309 550 2009 \n",
       "Q 550 2706 781 3412 \n",
       "Q 1013 4119 1484 4856 \n",
       "L 1984 4856 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-47\" d=\"M 3809 666 \n",
       "L 3809 1919 \n",
       "L 2778 1919 \n",
       "L 2778 2438 \n",
       "L 4434 2438 \n",
       "L 4434 434 \n",
       "Q 4069 175 3628 42 \n",
       "Q 3188 -91 2688 -91 \n",
       "Q 1594 -91 976 548 \n",
       "Q 359 1188 359 2328 \n",
       "Q 359 3472 976 4111 \n",
       "Q 1594 4750 2688 4750 \n",
       "Q 3144 4750 3555 4637 \n",
       "Q 3966 4525 4313 4306 \n",
       "L 4313 3634 \n",
       "Q 3963 3931 3569 4081 \n",
       "Q 3175 4231 2741 4231 \n",
       "Q 1884 4231 1454 3753 \n",
       "Q 1025 3275 1025 2328 \n",
       "Q 1025 1384 1454 906 \n",
       "Q 1884 428 2741 428 \n",
       "Q 3075 428 3337 486 \n",
       "Q 3600 544 3809 666 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-48\" d=\"M 628 4666 \n",
       "L 1259 4666 \n",
       "L 1259 2753 \n",
       "L 3553 2753 \n",
       "L 3553 4666 \n",
       "L 4184 4666 \n",
       "L 4184 0 \n",
       "L 3553 0 \n",
       "L 3553 2222 \n",
       "L 1259 2222 \n",
       "L 1259 0 \n",
       "L 628 0 \n",
       "L 628 4666 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-7a\" d=\"M 353 3500 \n",
       "L 3084 3500 \n",
       "L 3084 2975 \n",
       "L 922 459 \n",
       "L 3084 459 \n",
       "L 3084 0 \n",
       "L 275 0 \n",
       "L 275 525 \n",
       "L 2438 3041 \n",
       "L 353 3041 \n",
       "L 353 3500 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-29\" d=\"M 513 4856 \n",
       "L 1013 4856 \n",
       "Q 1481 4119 1714 3412 \n",
       "Q 1947 2706 1947 2009 \n",
       "Q 1947 1309 1714 600 \n",
       "Q 1481 -109 1013 -844 \n",
       "L 513 -844 \n",
       "Q 928 -128 1133 580 \n",
       "Q 1338 1288 1338 2009 \n",
       "Q 1338 2731 1133 3434 \n",
       "Q 928 4138 513 4856 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "      </defs>\n",
       "      <use xlink:href=\"#DejaVuSans-46\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-72\" x=\"50.269531\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-65\" x=\"89.132812\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-71\" x=\"150.65625\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-75\" x=\"214.132812\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-65\" x=\"277.511719\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-6e\" x=\"339.035156\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-63\" x=\"402.414062\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-79\" x=\"457.394531\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-20\" x=\"516.574219\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-28\" x=\"548.361328\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-47\" x=\"587.375\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-48\" x=\"664.865234\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-7a\" x=\"740.060547\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-29\" x=\"792.550781\"/>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "   <g id=\"matplotlib.axis_2\">\n",
       "    <g id=\"ytick_1\">\n",
       "     <g id=\"line2d_7\">\n",
       "      <defs>\n",
       "       <path id=\"m48afb21bd5\" d=\"M 0 0 \n",
       "L -3.5 0 \n",
       "\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </defs>\n",
       "      <g>\n",
       "       <use xlink:href=\"#m48afb21bd5\" x=\"58.471875\" y=\"319.638936\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_8\">\n",
       "      <!-- 0 -->\n",
       "      <g transform=\"translate(42.564375 324.957842)scale(0.14 -0.14)\">\n",
       "       <use xlink:href=\"#DejaVuSans-30\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_2\">\n",
       "     <g id=\"line2d_8\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#m48afb21bd5\" x=\"58.471875\" y=\"283.154009\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_9\">\n",
       "      <!-- 50 -->\n",
       "      <g transform=\"translate(33.656875 288.472915)scale(0.14 -0.14)\">\n",
       "       <use xlink:href=\"#DejaVuSans-35\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"63.623047\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_3\">\n",
       "     <g id=\"line2d_9\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#m48afb21bd5\" x=\"58.471875\" y=\"246.669082\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_10\">\n",
       "      <!-- 100 -->\n",
       "      <g transform=\"translate(24.749375 251.987988)scale(0.14 -0.14)\">\n",
       "       <use xlink:href=\"#DejaVuSans-31\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_4\">\n",
       "     <g id=\"line2d_10\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#m48afb21bd5\" x=\"58.471875\" y=\"210.184155\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_11\">\n",
       "      <!-- 150 -->\n",
       "      <g transform=\"translate(24.749375 215.503061)scale(0.14 -0.14)\">\n",
       "       <use xlink:href=\"#DejaVuSans-31\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-35\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_5\">\n",
       "     <g id=\"line2d_11\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#m48afb21bd5\" x=\"58.471875\" y=\"173.699228\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_12\">\n",
       "      <!-- 200 -->\n",
       "      <g transform=\"translate(24.749375 179.018135)scale(0.14 -0.14)\">\n",
       "       <use xlink:href=\"#DejaVuSans-32\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_6\">\n",
       "     <g id=\"line2d_12\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#m48afb21bd5\" x=\"58.471875\" y=\"137.214302\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_13\">\n",
       "      <!-- 250 -->\n",
       "      <g transform=\"translate(24.749375 142.533208)scale(0.14 -0.14)\">\n",
       "       <use xlink:href=\"#DejaVuSans-32\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-35\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_7\">\n",
       "     <g id=\"line2d_13\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#m48afb21bd5\" x=\"58.471875\" y=\"100.729375\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_14\">\n",
       "      <!-- 300 -->\n",
       "      <g transform=\"translate(24.749375 106.048281)scale(0.14 -0.14)\">\n",
       "       <use xlink:href=\"#DejaVuSans-33\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_8\">\n",
       "     <g id=\"line2d_14\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#m48afb21bd5\" x=\"58.471875\" y=\"64.244448\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_15\">\n",
       "      <!-- 350 -->\n",
       "      <g transform=\"translate(24.749375 69.563354)scale(0.14 -0.14)\">\n",
       "       <use xlink:href=\"#DejaVuSans-33\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-35\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"ytick_9\">\n",
       "     <g id=\"line2d_15\">\n",
       "      <g>\n",
       "       <use xlink:href=\"#m48afb21bd5\" x=\"58.471875\" y=\"27.759521\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "     <g id=\"text_16\">\n",
       "      <!-- 400 -->\n",
       "      <g transform=\"translate(24.749375 33.078428)scale(0.14 -0.14)\">\n",
       "       <use xlink:href=\"#DejaVuSans-34\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"63.623047\"/>\n",
       "       <use xlink:href=\"#DejaVuSans-30\" x=\"127.246094\"/>\n",
       "      </g>\n",
       "     </g>\n",
       "    </g>\n",
       "    <g id=\"text_17\">\n",
       "     <!-- Intensity (MJy/sr) -->\n",
       "     <g transform=\"translate(17.837813 237.457906)rotate(-90)scale(0.14 -0.14)\">\n",
       "      <defs>\n",
       "       <path id=\"DejaVuSans-49\" d=\"M 628 4666 \n",
       "L 1259 4666 \n",
       "L 1259 0 \n",
       "L 628 0 \n",
       "L 628 4666 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-74\" d=\"M 1172 4494 \n",
       "L 1172 3500 \n",
       "L 2356 3500 \n",
       "L 2356 3053 \n",
       "L 1172 3053 \n",
       "L 1172 1153 \n",
       "Q 1172 725 1289 603 \n",
       "Q 1406 481 1766 481 \n",
       "L 2356 481 \n",
       "L 2356 0 \n",
       "L 1766 0 \n",
       "Q 1100 0 847 248 \n",
       "Q 594 497 594 1153 \n",
       "L 594 3053 \n",
       "L 172 3053 \n",
       "L 172 3500 \n",
       "L 594 3500 \n",
       "L 594 4494 \n",
       "L 1172 4494 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-73\" d=\"M 2834 3397 \n",
       "L 2834 2853 \n",
       "Q 2591 2978 2328 3040 \n",
       "Q 2066 3103 1784 3103 \n",
       "Q 1356 3103 1142 2972 \n",
       "Q 928 2841 928 2578 \n",
       "Q 928 2378 1081 2264 \n",
       "Q 1234 2150 1697 2047 \n",
       "L 1894 2003 \n",
       "Q 2506 1872 2764 1633 \n",
       "Q 3022 1394 3022 966 \n",
       "Q 3022 478 2636 193 \n",
       "Q 2250 -91 1575 -91 \n",
       "Q 1294 -91 989 -36 \n",
       "Q 684 19 347 128 \n",
       "L 347 722 \n",
       "Q 666 556 975 473 \n",
       "Q 1284 391 1588 391 \n",
       "Q 1994 391 2212 530 \n",
       "Q 2431 669 2431 922 \n",
       "Q 2431 1156 2273 1281 \n",
       "Q 2116 1406 1581 1522 \n",
       "L 1381 1569 \n",
       "Q 847 1681 609 1914 \n",
       "Q 372 2147 372 2553 \n",
       "Q 372 3047 722 3315 \n",
       "Q 1072 3584 1716 3584 \n",
       "Q 2034 3584 2315 3537 \n",
       "Q 2597 3491 2834 3397 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-69\" d=\"M 603 3500 \n",
       "L 1178 3500 \n",
       "L 1178 0 \n",
       "L 603 0 \n",
       "L 603 3500 \n",
       "z\n",
       "M 603 4863 \n",
       "L 1178 4863 \n",
       "L 1178 4134 \n",
       "L 603 4134 \n",
       "L 603 4863 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-4d\" d=\"M 628 4666 \n",
       "L 1569 4666 \n",
       "L 2759 1491 \n",
       "L 3956 4666 \n",
       "L 4897 4666 \n",
       "L 4897 0 \n",
       "L 4281 0 \n",
       "L 4281 4097 \n",
       "L 3078 897 \n",
       "L 2444 897 \n",
       "L 1241 4097 \n",
       "L 1241 0 \n",
       "L 628 0 \n",
       "L 628 4666 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-4a\" d=\"M 628 4666 \n",
       "L 1259 4666 \n",
       "L 1259 325 \n",
       "Q 1259 -519 939 -900 \n",
       "Q 619 -1281 -91 -1281 \n",
       "L -331 -1281 \n",
       "L -331 -750 \n",
       "L -134 -750 \n",
       "Q 284 -750 456 -515 \n",
       "Q 628 -281 628 325 \n",
       "L 628 4666 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "       <path id=\"DejaVuSans-2f\" d=\"M 1625 4666 \n",
       "L 2156 4666 \n",
       "L 531 -594 \n",
       "L 0 -594 \n",
       "L 1625 4666 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "      </defs>\n",
       "      <use xlink:href=\"#DejaVuSans-49\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-6e\" x=\"29.492188\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-74\" x=\"92.871094\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-65\" x=\"132.080078\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-6e\" x=\"193.603516\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-73\" x=\"256.982422\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-69\" x=\"309.082031\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-74\" x=\"336.865234\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-79\" x=\"376.074219\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-20\" x=\"435.253906\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-28\" x=\"467.041016\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-4d\" x=\"506.054688\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-4a\" x=\"592.333984\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-79\" x=\"621.826172\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-2f\" x=\"681.005859\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-73\" x=\"714.697266\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-72\" x=\"766.796875\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-29\" x=\"807.910156\"/>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "   <g id=\"LineCollection_1\">\n",
       "    <path d=\"M 76.22642 173.181872 \n",
       "L 76.22642 173.161441 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 84.609994 137.587178 \n",
       "L 84.609994 137.559449 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 93.179869 105.837994 \n",
       "L 93.179869 105.801509 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 101.563442 80.48243 \n",
       "L 101.563442 80.448863 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 109.947016 61.282602 \n",
       "L 109.947016 61.250495 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 118.516891 48.148758 \n",
       "L 118.516891 48.11811 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 126.900464 41.277187 \n",
       "L 126.900464 41.250918 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 135.470339 39.828735 \n",
       "L 135.470339 39.802466 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 143.853913 43.167106 \n",
       "L 143.853913 43.143756 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 152.237486 50.512251 \n",
       "L 152.237486 50.49182 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 160.807361 61.289169 \n",
       "L 160.807361 61.270197 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 169.190934 74.266128 \n",
       "L 169.190934 74.248615 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 177.574508 89.006768 \n",
       "L 177.574508 88.990715 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 186.144383 105.17032 \n",
       "L 186.144383 105.155726 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 194.527956 121.583429 \n",
       "L 194.527956 121.567376 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 203.097831 138.508057 \n",
       "L 203.097831 138.490544 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 211.481405 154.781064 \n",
       "L 211.481405 154.760633 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 219.864978 170.553498 \n",
       "L 219.864978 170.530148 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 228.434853 185.926057 \n",
       "L 228.434853 185.899788 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 236.818427 200.108478 \n",
       "L 236.818427 200.076371 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 245.202 213.301427 \n",
       "L 245.202 213.269321 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 253.771875 225.645008 \n",
       "L 253.771875 225.611442 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 262.155448 236.785315 \n",
       "L 262.155448 236.751749 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 270.725323 247.08647 \n",
       "L 270.725323 247.052903 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 279.108897 256.144947 \n",
       "L 279.108897 256.11284 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 287.49247 264.287653 \n",
       "L 287.49247 264.257006 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 296.062345 271.664176 \n",
       "L 296.062345 271.634988 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 304.445919 278.054146 \n",
       "L 304.445919 278.026417 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 312.829492 283.734849 \n",
       "L 312.829492 283.70712 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 321.399367 288.810632 \n",
       "L 321.399367 288.782903 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 329.782941 293.128258 \n",
       "L 329.782941 293.097611 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 338.352816 296.989823 \n",
       "L 338.352816 296.956257 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 346.736389 300.262521 \n",
       "L 346.736389 300.224576 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 355.119962 303.136073 \n",
       "L 355.119962 303.09521 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 363.689837 305.610481 \n",
       "L 363.689837 305.566699 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 372.073411 307.701797 \n",
       "L 372.073411 307.655097 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 380.643286 309.585149 \n",
       "L 380.643286 309.536989 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 389.026859 311.115327 \n",
       "L 389.026859 311.064248 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 397.410433 312.429514 \n",
       "L 397.410433 312.369679 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 405.980308 313.57587 \n",
       "L 405.980308 313.495604 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 414.363881 314.531776 \n",
       "L 414.363881 314.403349 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 422.747455 315.519058 \n",
       "L 422.747455 315.292851 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "    <path d=\"M 431.31733 316.544284 \n",
       "L 431.31733 316.132734 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke: #000000; stroke-width: 1.5\"/>\n",
       "   </g>\n",
       "   <g id=\"line2d_16\">\n",
       "    <path d=\"M 76.22642 173.176891 \n",
       "L 84.609994 137.582029 \n",
       "L 93.179869 105.833427 \n",
       "L 101.563442 80.471707 \n",
       "L 109.947016 61.284084 \n",
       "L 118.516891 48.115983 \n",
       "L 126.900464 41.246639 \n",
       "L 135.470339 39.813449 \n",
       "L 143.853913 43.183663 \n",
       "L 152.237486 50.510498 \n",
       "L 160.807361 61.283349 \n",
       "L 169.190934 74.272183 \n",
       "L 177.574508 88.98782 \n",
       "L 186.144383 105.174053 \n",
       "L 194.527956 121.586145 \n",
       "L 203.097831 138.488555 \n",
       "L 211.481405 154.7818 \n",
       "L 219.864978 170.554967 \n",
       "L 228.434853 185.924946 \n",
       "L 236.818427 200.075013 \n",
       "L 245.202 213.255071 \n",
       "L 253.771875 225.673685 \n",
       "L 262.155448 236.775743 \n",
       "L 270.725323 247.067667 \n",
       "L 279.108897 256.135524 \n",
       "L 287.49247 264.262017 \n",
       "L 296.062345 271.65609 \n",
       "L 304.445919 278.061062 \n",
       "L 312.829492 283.713827 \n",
       "L 321.399367 288.784605 \n",
       "L 329.782941 293.119973 \n",
       "L 338.352816 296.978928 \n",
       "L 346.736389 300.254515 \n",
       "L 355.119962 303.092181 \n",
       "L 363.689837 305.593647 \n",
       "L 372.073411 307.697789 \n",
       "L 380.643286 309.542478 \n",
       "L 389.026859 311.086163 \n",
       "L 397.410433 312.405919 \n",
       "L 405.980308 313.554799 \n",
       "L 414.363881 314.509816 \n",
       "L 422.747455 315.321217 \n",
       "L 431.31733 316.023353 \n",
       "\" clip-path=\"url(#pa6b7f98970)\" style=\"fill: none; stroke-dasharray: 5.55,2.4; stroke-dashoffset: 0; stroke: #06c2ac; stroke-width: 1.5\"/>\n",
       "   </g>\n",
       "   <g id=\"line2d_17\">\n",
       "    <defs>\n",
       "     <path id=\"mf8a02c760c\" d=\"M 0 1.5 \n",
       "C 0.397805 1.5 0.77937 1.341951 1.06066 1.06066 \n",
       "C 1.341951 0.77937 1.5 0.397805 1.5 0 \n",
       "C 1.5 -0.397805 1.341951 -0.77937 1.06066 -1.06066 \n",
       "C 0.77937 -1.341951 0.397805 -1.5 0 -1.5 \n",
       "C -0.397805 -1.5 -0.77937 -1.341951 -1.06066 -1.06066 \n",
       "C -1.341951 -0.77937 -1.5 -0.397805 -1.5 0 \n",
       "C -1.5 0.397805 -1.341951 0.77937 -1.06066 1.06066 \n",
       "C -0.77937 1.341951 -0.397805 1.5 0 1.5 \n",
       "z\n",
       "\" style=\"stroke: #000000\"/>\n",
       "    </defs>\n",
       "    <g clip-path=\"url(#pa6b7f98970)\">\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"76.22642\" y=\"173.171656\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"84.609994\" y=\"137.573313\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"93.179869\" y=\"105.819752\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"101.563442\" y=\"80.465647\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"109.947016\" y=\"61.266548\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"118.516891\" y=\"48.133434\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"126.900464\" y=\"41.264052\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"135.470339\" y=\"39.8156\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"143.853913\" y=\"43.155431\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"152.237486\" y=\"50.502036\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"160.807361\" y=\"61.279683\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"169.190934\" y=\"74.257371\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"177.574508\" y=\"88.998741\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"186.144383\" y=\"105.163023\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"194.527956\" y=\"121.575403\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"203.097831\" y=\"138.499301\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"211.481405\" y=\"154.770848\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"219.864978\" y=\"170.541823\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"228.434853\" y=\"185.912923\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"236.818427\" y=\"200.092424\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"245.202\" y=\"213.285374\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"253.771875\" y=\"225.628225\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"262.155448\" y=\"236.768532\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"270.725323\" y=\"247.069687\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"279.108897\" y=\"256.128894\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"287.49247\" y=\"264.272329\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"296.062345\" y=\"271.649582\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"304.445919\" y=\"278.040281\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"312.829492\" y=\"283.720985\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"321.399367\" y=\"288.796768\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"329.782941\" y=\"293.112934\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"338.352816\" y=\"296.97304\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"346.736389\" y=\"300.243549\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"355.119962\" y=\"303.115642\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"363.689837\" y=\"305.58859\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"372.073411\" y=\"307.678447\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"380.643286\" y=\"309.561069\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"389.026859\" y=\"311.089788\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"397.410433\" y=\"312.399596\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"405.980308\" y=\"313.535737\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"414.363881\" y=\"314.467562\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"422.747455\" y=\"315.405954\" style=\"stroke: #000000\"/>\n",
       "     <use xlink:href=\"#mf8a02c760c\" x=\"431.31733\" y=\"316.338509\" style=\"stroke: #000000\"/>\n",
       "    </g>\n",
       "   </g>\n",
       "   <g id=\"patch_3\">\n",
       "    <path d=\"M 58.471875 330.381375 \n",
       "L 58.471875 25.965375 \n",
       "\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_4\">\n",
       "    <path d=\"M 449.071875 330.381375 \n",
       "L 449.071875 25.965375 \n",
       "\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_5\">\n",
       "    <path d=\"M 58.471875 330.381375 \n",
       "L 449.071875 330.381375 \n",
       "\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_6\">\n",
       "    <path d=\"M 58.471875 25.965375 \n",
       "L 449.071875 25.965375 \n",
       "\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n",
       "   </g>\n",
       "   <g id=\"text_18\">\n",
       "    <!-- Monopole temperature -->\n",
       "    <g transform=\"translate(157.481625 19.965375)scale(0.168 -0.168)\">\n",
       "     <defs>\n",
       "      <path id=\"DejaVuSans-6f\" d=\"M 1959 3097 \n",
       "Q 1497 3097 1228 2736 \n",
       "Q 959 2375 959 1747 \n",
       "Q 959 1119 1226 758 \n",
       "Q 1494 397 1959 397 \n",
       "Q 2419 397 2687 759 \n",
       "Q 2956 1122 2956 1747 \n",
       "Q 2956 2369 2687 2733 \n",
       "Q 2419 3097 1959 3097 \n",
       "z\n",
       "M 1959 3584 \n",
       "Q 2709 3584 3137 3096 \n",
       "Q 3566 2609 3566 1747 \n",
       "Q 3566 888 3137 398 \n",
       "Q 2709 -91 1959 -91 \n",
       "Q 1206 -91 779 398 \n",
       "Q 353 888 353 1747 \n",
       "Q 353 2609 779 3096 \n",
       "Q 1206 3584 1959 3584 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "      <path id=\"DejaVuSans-70\" d=\"M 1159 525 \n",
       "L 1159 -1331 \n",
       "L 581 -1331 \n",
       "L 581 3500 \n",
       "L 1159 3500 \n",
       "L 1159 2969 \n",
       "Q 1341 3281 1617 3432 \n",
       "Q 1894 3584 2278 3584 \n",
       "Q 2916 3584 3314 3078 \n",
       "Q 3713 2572 3713 1747 \n",
       "Q 3713 922 3314 415 \n",
       "Q 2916 -91 2278 -91 \n",
       "Q 1894 -91 1617 61 \n",
       "Q 1341 213 1159 525 \n",
       "z\n",
       "M 3116 1747 \n",
       "Q 3116 2381 2855 2742 \n",
       "Q 2594 3103 2138 3103 \n",
       "Q 1681 3103 1420 2742 \n",
       "Q 1159 2381 1159 1747 \n",
       "Q 1159 1113 1420 752 \n",
       "Q 1681 391 2138 391 \n",
       "Q 2594 391 2855 752 \n",
       "Q 3116 1113 3116 1747 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "      <path id=\"DejaVuSans-6c\" d=\"M 603 4863 \n",
       "L 1178 4863 \n",
       "L 1178 0 \n",
       "L 603 0 \n",
       "L 603 4863 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "      <path id=\"DejaVuSans-6d\" d=\"M 3328 2828 \n",
       "Q 3544 3216 3844 3400 \n",
       "Q 4144 3584 4550 3584 \n",
       "Q 5097 3584 5394 3201 \n",
       "Q 5691 2819 5691 2113 \n",
       "L 5691 0 \n",
       "L 5113 0 \n",
       "L 5113 2094 \n",
       "Q 5113 2597 4934 2840 \n",
       "Q 4756 3084 4391 3084 \n",
       "Q 3944 3084 3684 2787 \n",
       "Q 3425 2491 3425 1978 \n",
       "L 3425 0 \n",
       "L 2847 0 \n",
       "L 2847 2094 \n",
       "Q 2847 2600 2669 2842 \n",
       "Q 2491 3084 2119 3084 \n",
       "Q 1678 3084 1418 2786 \n",
       "Q 1159 2488 1159 1978 \n",
       "L 1159 0 \n",
       "L 581 0 \n",
       "L 581 3500 \n",
       "L 1159 3500 \n",
       "L 1159 2956 \n",
       "Q 1356 3278 1631 3431 \n",
       "Q 1906 3584 2284 3584 \n",
       "Q 2666 3584 2933 3390 \n",
       "Q 3200 3197 3328 2828 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "      <path id=\"DejaVuSans-61\" d=\"M 2194 1759 \n",
       "Q 1497 1759 1228 1600 \n",
       "Q 959 1441 959 1056 \n",
       "Q 959 750 1161 570 \n",
       "Q 1363 391 1709 391 \n",
       "Q 2188 391 2477 730 \n",
       "Q 2766 1069 2766 1631 \n",
       "L 2766 1759 \n",
       "L 2194 1759 \n",
       "z\n",
       "M 3341 1997 \n",
       "L 3341 0 \n",
       "L 2766 0 \n",
       "L 2766 531 \n",
       "Q 2569 213 2275 61 \n",
       "Q 1981 -91 1556 -91 \n",
       "Q 1019 -91 701 211 \n",
       "Q 384 513 384 1019 \n",
       "Q 384 1609 779 1909 \n",
       "Q 1175 2209 1959 2209 \n",
       "L 2766 2209 \n",
       "L 2766 2266 \n",
       "Q 2766 2663 2505 2880 \n",
       "Q 2244 3097 1772 3097 \n",
       "Q 1472 3097 1187 3025 \n",
       "Q 903 2953 641 2809 \n",
       "L 641 3341 \n",
       "Q 956 3463 1253 3523 \n",
       "Q 1550 3584 1831 3584 \n",
       "Q 2591 3584 2966 3190 \n",
       "Q 3341 2797 3341 1997 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "     </defs>\n",
       "     <use xlink:href=\"#DejaVuSans-4d\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-6f\" x=\"86.279297\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-6e\" x=\"147.460938\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-6f\" x=\"210.839844\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-70\" x=\"272.021484\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-6f\" x=\"335.498047\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-6c\" x=\"396.679688\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-65\" x=\"424.462891\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-20\" x=\"485.986328\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-74\" x=\"517.773438\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-65\" x=\"556.982422\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-6d\" x=\"618.505859\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-70\" x=\"715.917969\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-65\" x=\"779.394531\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-72\" x=\"840.917969\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-61\" x=\"882.03125\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-74\" x=\"943.310547\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-75\" x=\"982.519531\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-72\" x=\"1045.898438\"/>\n",
       "     <use xlink:href=\"#DejaVuSans-65\" x=\"1084.761719\"/>\n",
       "    </g>\n",
       "   </g>\n",
       "   <g id=\"legend_1\">\n",
       "    <g id=\"line2d_18\">\n",
       "     <path d=\"M 382.965625 44.303187 \n",
       "L 396.965625 44.303187 \n",
       "L 410.965625 44.303187 \n",
       "\" style=\"fill: none; stroke-dasharray: 5.55,2.4; stroke-dashoffset: 0; stroke: #06c2ac; stroke-width: 1.5\"/>\n",
       "    </g>\n",
       "    <g id=\"text_19\">\n",
       "     <!-- fit -->\n",
       "     <g transform=\"translate(422.165625 49.203187)scale(0.14 -0.14)\">\n",
       "      <defs>\n",
       "       <path id=\"DejaVuSans-66\" d=\"M 2375 4863 \n",
       "L 2375 4384 \n",
       "L 1825 4384 \n",
       "Q 1516 4384 1395 4259 \n",
       "Q 1275 4134 1275 3809 \n",
       "L 1275 3500 \n",
       "L 2222 3500 \n",
       "L 2222 3053 \n",
       "L 1275 3053 \n",
       "L 1275 0 \n",
       "L 697 0 \n",
       "L 697 3053 \n",
       "L 147 3053 \n",
       "L 147 3500 \n",
       "L 697 3500 \n",
       "L 697 3744 \n",
       "Q 697 4328 969 4595 \n",
       "Q 1241 4863 1831 4863 \n",
       "L 2375 4863 \n",
       "z\n",
       "\" transform=\"scale(0.015625)\"/>\n",
       "      </defs>\n",
       "      <use xlink:href=\"#DejaVuSans-66\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-69\" x=\"35.205078\"/>\n",
       "      <use xlink:href=\"#DejaVuSans-74\" x=\"62.988281\"/>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "  </g>\n",
       " </g>\n",
       " <defs>\n",
       "  <clipPath id=\"pa6b7f98970\">\n",
       "   <rect x=\"58.471875\" y=\"25.965375\" width=\"390.6\" height=\"304.416\"/>\n",
       "  </clipPath>\n",
       " </defs>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<Figure size 504x403.2 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bb_simple = BlackBody(2.725*u.K) #from the temperature given in the data description\n",
    "plt.figure()\n",
    "plt.errorbar(firas_freq, firas[\"monopole\"], yerr = firas[\"sigma\"], fmt = 'k.')\n",
    "plt.plot(firas_freq.to(\"GHz\"), bb_simple(firas_freq).to('MJy/sr'), color = 'xkcd:turquoise', ls = \"--\", label = \"fit\")\n",
    "plt.xlabel(\"Frequency (GHz)\")\n",
    "plt.ylabel(\"Intensity (MJy/sr)\")\n",
    "plt.title(\"Monopole temperature\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c7bff0c",
   "metadata": {},
   "source": [
    "### Consistency check 1\n",
    "\n",
    "The error bars are smaller than the fit line.\n",
    "\n",
    "1. Plot the error bars scaled by a sufficiently large factor.\n",
    "2. Note the frequencies where the error bars are noticeably larger. What could cause this relative increase?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dc422be",
   "metadata": {},
   "source": [
    "## 2. Fitting the data\n",
    "\n",
    "How do we know that this temperature is actually the most likely value? We can use `emcee` to reproduce the temperature used to generate the monopole data. If you've never used `emcee` before, a tutorial can be found in the [documentation](https://emcee.readthedocs.io/en/stable/tutorials/line/)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9924bb57-6db8-4135-b92a-a48e3974c578",
   "metadata": {},
   "source": [
    "We will first define the log-likelihood function for our model assessment. This will be passed into `emcee`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93b98098",
   "metadata": {},
   "source": [
    "### Consistency check 2\n",
    "\n",
    "The blackbody spectrum is given mathematically by\n",
    "$$\n",
    "S_\\nu(T)=\\frac{2h\\nu^3}{c^2}\\frac1{e^{h\\nu/kT}-1}.\n",
    "$$\n",
    "Assuming that the noise $S_\\mathrm{err}$ describes Gaussian noise, such that\n",
    "$$\n",
    "d_\\nu=S_\\nu(T) + n_\\nu\n",
    "$$\n",
    "where $n_\\nu\\sim\\mathcal N(0,S_\\mathrm{err}^2)$, derive the log-likelihood and code it in the next section.\n",
    "\n",
    "\n",
    "Potential pitfalls:\n",
    "\n",
    "1. Make sure that $h\\nu/kT$ is unitless.\n",
    "2. Make sure that the model and the data have the same units, MJy/sr.\n",
    "3. $e^{h\\nu/kT}-1$ can be very close to zero. The numpy function `np.expm1` may be useful for numerical stability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "84069fc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#log likelihood function:\n",
    "def log_LBB(T, nu, S, Serr):\n",
    "    #T is temperature in Kelvin\n",
    "    #nu is frequency in s-1\n",
    "    #S is the monopole spectrum in MJy/sr\n",
    "    #and Serr is the error in kJy/sr\n",
    "    \n",
    "    return log_L"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be15320c-580d-4941-b791-898cb0b9aece",
   "metadata": {},
   "source": [
    "Then we will run the `emcee` chain! We will need to define a few things to get `emcee` to run:\n",
    "A random seed, the chain length, as well as a initial position for the sampler to explore from, as well as "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9b74b914",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "emcee: Exception while calling your likelihood function:\n",
      "  params: [2.70049671] K\n",
      "  args: (<Quantity [ 68.05288797,  81.54354858,  95.33400164, 108.82466225,\n",
      "           122.31532286, 136.10577593, 149.59643654, 163.38688961,\n",
      "           176.87755022, 190.36821083, 204.1586639 , 217.64932451,\n",
      "           231.13998512, 244.93043819, 258.4210988 , 272.21155186,\n",
      "           285.70221247, 299.19287308, 312.98332615, 326.47398676,\n",
      "           339.96464737, 353.75510044, 367.24576105, 381.03621412,\n",
      "           394.52687473, 408.01753534, 421.80798841, 435.29864902,\n",
      "           448.78930963, 462.57976269, 476.0704233 , 489.86087637,\n",
      "           503.35153698, 516.84219759, 530.63265066, 544.12331127,\n",
      "           557.91376434, 571.40442495, 584.89508556, 598.68553863,\n",
      "           612.17619924, 625.66685985, 639.45731291] GHz>, <Quantity [200.723, 249.508, 293.024, 327.77 , 354.081, 372.079, 381.493,\n",
      "           383.478, 378.901, 368.833, 354.063, 336.278, 316.076, 293.924,\n",
      "           271.432, 248.239, 225.94 , 204.327, 183.262, 163.83 , 145.75 ,\n",
      "           128.835, 113.568,  99.451,  87.036,  75.876,  65.766,  57.008,\n",
      "            49.223,  42.267,  36.352,  31.062,  26.58 ,  22.644,  19.255,\n",
      "            16.391,  13.811,  11.716,   9.921,   8.364,   7.087,   5.801,\n",
      "             4.523] MJy / sr>, <Quantity [ 14.,  19.,  25.,  23.,  22.,  21.,  18.,  18.,  16.,  14.,\n",
      "            13.,  12.,  11.,  10.,  11.,  12.,  14.,  16.,  18.,  22.,\n",
      "            22.,  23.,  23.,  23.,  22.,  21.,  20.,  19.,  19.,  19.,\n",
      "            21.,  23.,  26.,  28.,  30.,  32.,  33.,  35.,  41.,  55.,\n",
      "            88., 155., 282.] kJy / sr>)\n",
      "  kwargs: {}\n",
      "  exception:\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Traceback (most recent call last):\n",
      "  File \"/Users/Katie/opt/anaconda3/lib/python3.9/site-packages/emcee/ensemble.py\", line 624, in __call__\n",
      "    return self.f(x, *self.args, **self.kwargs)\n",
      "  File \"/var/folders/cy/62nzdxhj0bs44z09p5_k64bm0000gn/T/ipykernel_64400/2755273366.py\", line 8, in log_LBB\n",
      "    return log_L\n",
      "NameError: name 'log_L' is not defined\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'log_L' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Input \u001b[0;32mIn [5]\u001b[0m, in \u001b[0;36m<cell line: 11>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[38;5;66;03m# Run the sampler\u001b[39;00m\n\u001b[1;32m     10\u001b[0m sampler \u001b[38;5;241m=\u001b[39m mc\u001b[38;5;241m.\u001b[39mEnsembleSampler(nwalkers, ndim, log_LBB, args \u001b[38;5;241m=\u001b[39m (firas_freq, firas[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmonopole\u001b[39m\u001b[38;5;124m\"\u001b[39m], firas[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124msigma\u001b[39m\u001b[38;5;124m\"\u001b[39m]))\n\u001b[0;32m---> 11\u001b[0m \u001b[43msampler\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun_mcmc\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpos\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mchain_len\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mprogress\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/emcee/ensemble.py:443\u001b[0m, in \u001b[0;36mEnsembleSampler.run_mcmc\u001b[0;34m(self, initial_state, nsteps, **kwargs)\u001b[0m\n\u001b[1;32m    440\u001b[0m     initial_state \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_previous_state\n\u001b[1;32m    442\u001b[0m results \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[0;32m--> 443\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m results \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msample(initial_state, iterations\u001b[38;5;241m=\u001b[39mnsteps, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[1;32m    444\u001b[0m     \u001b[38;5;28;01mpass\u001b[39;00m\n\u001b[1;32m    446\u001b[0m \u001b[38;5;66;03m# Store so that the ``initial_state=None`` case will work\u001b[39;00m\n",
      "File \u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/emcee/ensemble.py:344\u001b[0m, in \u001b[0;36mEnsembleSampler.sample\u001b[0;34m(self, initial_state, log_prob0, rstate0, blobs0, iterations, tune, skip_initial_state_check, thin_by, thin, store, progress, progress_kwargs)\u001b[0m\n\u001b[1;32m    342\u001b[0m     state\u001b[38;5;241m.\u001b[39mblobs \u001b[38;5;241m=\u001b[39m blobs0\n\u001b[1;32m    343\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m state\u001b[38;5;241m.\u001b[39mlog_prob \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m--> 344\u001b[0m     state\u001b[38;5;241m.\u001b[39mlog_prob, state\u001b[38;5;241m.\u001b[39mblobs \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcompute_log_prob\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstate\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcoords\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    345\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m np\u001b[38;5;241m.\u001b[39mshape(state\u001b[38;5;241m.\u001b[39mlog_prob) \u001b[38;5;241m!=\u001b[39m (\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mnwalkers,):\n\u001b[1;32m    346\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mincompatible input dimensions\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "File \u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/emcee/ensemble.py:489\u001b[0m, in \u001b[0;36mEnsembleSampler.compute_log_prob\u001b[0;34m(self, coords)\u001b[0m\n\u001b[1;32m    487\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m    488\u001b[0m         map_func \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mmap\u001b[39m\n\u001b[0;32m--> 489\u001b[0m     results \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mlist\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mmap_func\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlog_prob_fn\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mp\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    491\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m    492\u001b[0m     log_prob \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray([\u001b[38;5;28mfloat\u001b[39m(l[\u001b[38;5;241m0\u001b[39m]) \u001b[38;5;28;01mfor\u001b[39;00m l \u001b[38;5;129;01min\u001b[39;00m results])\n",
      "File \u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/emcee/ensemble.py:624\u001b[0m, in \u001b[0;36m_FunctionWrapper.__call__\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m    622\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__call__\u001b[39m(\u001b[38;5;28mself\u001b[39m, x):\n\u001b[1;32m    623\u001b[0m     \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m--> 624\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mf\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    625\u001b[0m     \u001b[38;5;28;01mexcept\u001b[39;00m:  \u001b[38;5;66;03m# pragma: no cover\u001b[39;00m\n\u001b[1;32m    626\u001b[0m         \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mtraceback\u001b[39;00m\n",
      "Input \u001b[0;32mIn [4]\u001b[0m, in \u001b[0;36mlog_LBB\u001b[0;34m(T, nu, S, Serr)\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mlog_LBB\u001b[39m(T, nu, S, Serr):\n\u001b[1;32m      3\u001b[0m     \u001b[38;5;66;03m#T is temperature in Kelvin\u001b[39;00m\n\u001b[1;32m      4\u001b[0m     \u001b[38;5;66;03m#nu is frequency in s-1\u001b[39;00m\n\u001b[1;32m      5\u001b[0m     \u001b[38;5;66;03m#S is the monopole spectrum in MJy/sr\u001b[39;00m\n\u001b[1;32m      6\u001b[0m     \u001b[38;5;66;03m#and Serr is the error in kJy/sr\u001b[39;00m\n\u001b[0;32m----> 8\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mlog_L\u001b[49m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'log_L' is not defined"
     ]
    }
   ],
   "source": [
    "#run mcmc\n",
    "np.random.seed(42)\n",
    "chain_len = 3000\n",
    "T_ini = 2.7*u.K\n",
    "# Draw initial positions for each of the walkers\n",
    "pos = T_ini + 1e-3*np.random.randn(32, 1)*u.K\n",
    "nwalkers, ndim = pos.shape\n",
    "\n",
    "# Run the sampler\n",
    "sampler = mc.EnsembleSampler(nwalkers, ndim, log_LBB, args = (firas_freq, firas[\"monopole\"], firas[\"sigma\"]))\n",
    "sampler.run_mcmc(pos, chain_len, progress = True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f82baea",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get all samples:\n",
    "samples = sampler.get_chain()\n",
    "\n",
    "#plot samples over time\n",
    "plt.figure()\n",
    "plt.plot(samples[:, :, 0], \"k\", alpha = 0.3)\n",
    "plt.ylabel(\"T\")\n",
    "plt.xlabel(\"sample number\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c32d7a9",
   "metadata": {},
   "source": [
    "## Consistency check 3\n",
    "\n",
    "The autocorrelation length (see https://emcee.readthedocs.io/en/stable/tutorials/autocorr/) is an estimate for how many samples must be separated by to be considered independent. Using the estimate of the autocorrelation length $\\tau$, see how many integer multiples of $\\tau$ are necessary to see if the chain is truly burned in.\n",
    "\n",
    "If the chain is truly converged, a good rule of thumb is that it should look like a \"hairy caterpillar\". Comment on whether this seems like a good metric, and if you agree with it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5e2a302",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get autocorrelation time:\n",
    "tau = sampler.get_autocorr_time()\n",
    "print(tau)\n",
    "\n",
    "\n",
    "#plot samples over time\n",
    "plt.figure()\n",
    "plt.title('Burned in chain')\n",
    "plt.ylabel(\"T\")\n",
    "plt.xlabel(\"sample number\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab275d6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#print final estimate\n",
    "flat_samples = sampler.get_chain(discard = , thin = int(tau), flat = True)\n",
    "T_mcmc = np.percentile(flat_samples[:, 0], [16, 50, 84])\n",
    "q = np.diff(T_mcmc)\n",
    "\n",
    "txt = \"\\mathrm{{{3}}} = {0:.6f}_{{-{1:.6f}}}^{{{2:.6f}}}\"\n",
    "txt = txt.format(T_mcmc[1], q[0], q[1], \"T\")\n",
    "display(Math(txt))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc980fb2",
   "metadata": {},
   "source": [
    "## Consistency check 4\n",
    "\n",
    "Fixsen (2009) (https://ui.adsabs.harvard.edu/abs/2009ApJ...707..916F/abstract) is the widely-cited paper that contains the final temperature estimate of the FIRAS instrument. How does your result compare with these results? Comment on any notable differences, and posit a possible explanation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "050849b0",
   "metadata": {},
   "source": [
    "#### Short answer\n",
    "\n",
    "\n",
    "*I believe that...*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7adff4f9",
   "metadata": {},
   "source": [
    "# Spectral Distortions\n",
    "\n",
    "Now, we'll use the residuals to sample the spectral distortions, or the differences between the FIRAS data and a pure blackbody spectrum. These spectral distortions originate from energy that has been injected into the early universe, causing the system to be at a slightly higher energy state than that of pure equilibrium. We use two different parameters to quantify these distortions: $\\mu$ and $y$. These represent energy injections from high-energy scattering events at higher redshifts, and Compton scattering at lower redshifts, respectively (Fixsen et. al. 1996). \n",
    "\n",
    "Sunyaev and Zeldovich (1970) state that the source of $\\mu$-type spectral distortions is that prior to recombination and decoupling, at redshifts between $10^5 < z < 3\\times 10^6$, \"bremsstrahlung processes are not able to maintain complete thermodynamic equilibrium between matter and radiation. Meanwhile the transfer of energy between electrons and photons due to scattering occurs sufficiently quickly for a Bose-Einstein distribution to be established.\" So basically, electrons and such are being excited at a high enough rate that the energy bound up in them can't be radiated away fast enough for the plasma to reach equilibrium, creating a chemical potential $\\mu$.\n",
    "\n",
    "According to Fixsen et. al. (1996), the effect of $y$-distortions is to produce \"a Comptonized spectrum, [i.e.] a mixture of blackbodies at a range of temperatures.\" These typically come into play at $z<10^5$, and Chluba (2005) implies that they can also be created after decoupling by the Sunyaev-Zeldovich effect."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01dae0a9",
   "metadata": {},
   "source": [
    "### Consistency check 4\n",
    "\n",
    "At this point in the analysis, it is helpful to assume that we basically have the correct temperature, $T_0$, and any corrections will be small, and given by $\\Delta T$. This means that the data can be modeled as\n",
    "$$\n",
    "B_\\nu(T_0+\\Delta T)-B_\\nu(T_0)\\simeq\\,\\Delta T\n",
    "\\left.\\frac{\\partial B_\\nu}{\\partial T}\\right|_{T=T_0}\n",
    "$$\n",
    "\n",
    "Derive the spectrum of the deviations from the ideal blackbody, $\\partial B_\\nu/\\partial T$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "844507f3",
   "metadata": {},
   "source": [
    "*First, we...*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e71fa972-edd4-4f10-98f0-65a056210375",
   "metadata": {},
   "source": [
    "Now we want to sample for any potential distortions in the spectrum, whether $\\mu$ or $y$ type. Again we'll need to define log-likelihood functions for the spectral distortion distributions. Before we do that, let us define functional forms for both of the distortion types, $\\mu$ and $y$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfbbe62f",
   "metadata": {},
   "source": [
    "### Consistency check 5\n",
    "\n",
    "A $\\mu$-type distribution is equivalent to a Bose-Einstein distribution with chemical potential $\\mu$,\n",
    "$$\n",
    "I_\\nu^\\mathrm{BE}=\\frac{2h\\nu^3}{c^2}\\frac{1}{e^{x+\\mu}-1}\n",
    "$$\n",
    "where $x=h\\nu/kT$.\n",
    "\n",
    "Derive the first order correction, $I_\\nu^\\mathrm{SD,\\mu}=\\mu\\frac{\\partial I_\\nu}{\\partial \\mu}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d2a6341",
   "metadata": {},
   "source": [
    "*We begin by...*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "840bc1da",
   "metadata": {},
   "source": [
    " \n",
    " A similar approach can be followed to separate $y$ distortion map, where the emission law is\n",
    " \n",
    " \n",
    "\\begin{equation}\n",
    "I_{\\nu}^{\\mathrm{SD},y} = yT_0 \\Big(x\\frac{\\exp{x}+1}{\\exp{x}-1}-4\\Big)\\frac{\\partial B_\\nu}{\\partial T}\\Big|_{T=T_0}.\"\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08e6f206-90b8-4523-9746-67b679a2c75a",
   "metadata": {},
   "source": [
    "The above equations are the functions which we need to evaluate in order to fit for the spectral parameter. Let's do that now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5bfcd88",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dBdT(nu, T = 2.725*u.K):\n",
    "    #nu is frequency in s-1\n",
    "    x = const.h*nu/const.k_B/T\n",
    "    dBdT = \n",
    "    return dBdT\n",
    "\n",
    "\n",
    "def func_mu(nu, mu, T = 2.725):\n",
    "    #nu is frequency in s-1\n",
    "    #mu is dimensionless\n",
    "    #T is unitless or else emcee freaks out\n",
    "    T = T*u.K\n",
    "    x = const.h*nu/const.k_B/T\n",
    "    f = \n",
    "    return f\n",
    "    \n",
    "\n",
    "def func_y(nu, y, T = 2.725):\n",
    "    #nu is frequency in s-1\n",
    "    #y is dimensionless\n",
    "    #T is unitless or else emcee freaks out\n",
    "    T = T*u.K\n",
    "    x = const.h*nu/const.k_B/T\n",
    "    f = \n",
    "    return f"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7de3db5-382c-4b4f-905d-9b7a04644327",
   "metadata": {},
   "source": [
    "Define some likelihood functions here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6b2a104",
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_Ly(theta, nu, Sres, Serr):\n",
    "    #nu is an array of frequencies in s-1\n",
    "    #Sres and Serr are the residual and error arrays in kJy/sr\n",
    "    y = theta\n",
    "    model = func_y(nu, y)\n",
    "    \n",
    "    log_L = \n",
    "    return log_L\n",
    "\n",
    "def log_Lmu(theta, nu, Sres, Serr):\n",
    "    mu = theta\n",
    "    model = func_mu(nu, mu)\n",
    "    \n",
    "    log_L = \n",
    "    return log_L"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80967fb7-f680-46bc-abb2-9e972471f2e6",
   "metadata": {},
   "source": [
    "Before we try sampling for each of the parameters ($T$, $\\mu$, $y$) together, let's do them one at a time. Now setting up `emcee` for $y$: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96844fbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#mcmc for y:\n",
    "np.random.seed(42)\n",
    "chain_len = 2000\n",
    "y_ini = -1e-6\n",
    "pos_y = y_ini + 1e-6*np.random.randn(32, 1)\n",
    "nwalkers_y, ndim_y = pos_y.shape\n",
    "\n",
    "sampler_y = mc.EnsembleSampler(nwalkers_y, ndim_y, log_Ly, args = (firas_freq, firas[\"res\"], firas[\"sigma\"]))\n",
    "sampler_y.run_mcmc(pos_y, chain_len, progress = True);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a3fb611-4ecd-4589-a901-85f43eb5b3fe",
   "metadata": {},
   "source": [
    "Now for $\\mu$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c42c0b65",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_ini = -1e-5\n",
    "pos_mu = mu_ini + 1e-6*np.random.randn(32, 1)\n",
    "nwalkers_mu, ndim_mu = pos_mu.shape\n",
    "\n",
    "sampler_mu = mc.EnsembleSampler(nwalkers_mu, ndim_mu, log_Lmu, args = (firas_freq, firas[\"res\"], firas[\"sigma\"]))\n",
    "sampler_mu.run_mcmc(pos_mu, chain_len, progress = True);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea8de032-1cbc-45bf-b5c0-b06795d3bdd1",
   "metadata": {},
   "source": [
    "Let's plot the chains to visualize convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a41fa27c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#get autocorrelation time:\n",
    "tau_y = \n",
    "tau_mu = \n",
    "print(tau_y, tau_mu)\n",
    "\n",
    "#get all samples:\n",
    "samples_y = \n",
    "samples_mu = \n",
    "#print(samples_y.shape, samples_mu.shape)\n",
    "\n",
    "#plot samples over time\n",
    "fig, axes = plt.subplots(sharex=True, nrows=2)\n",
    "axes[0].plot()\n",
    "axes[0].set_ylabel(r\"$y$\")\n",
    "axes[1].plot()\n",
    "axes[1].set_ylabel(r\"$\\mu$\")\n",
    "axes[1].set_xlabel(\"sample number\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0bb1a2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "flat_samples_y = sampler_y.get_chain(discard = , thin = int(tau_y), flat = True)\n",
    "y_mcmc = np.percentile(flat_samples_y[:, 0], [16, 50, 84])\n",
    "q_y = np.diff(y_mcmc)\n",
    "\n",
    "flat_samples_mu = sampler_mu.get_chain(discard = , thin = int(tau_mu), flat = True)\n",
    "mu_mcmc = np.percentile(flat_samples_mu[:, 0], [16, 50, 84])\n",
    "q_mu = np.diff(mu_mcmc)\n",
    "\n",
    "txt = \"\\mathrm{{{3}}} = {0:.7f}_{{-{1:.6f}}}^{{{2:.6f}}}\"\n",
    "txt_y = txt.format(y_mcmc[1], q_y[0], q_y[1], \"y\")\n",
    "txt_mu = txt.format(mu_mcmc[1], q_mu[0], q_y[1], \"\\mu\")\n",
    "display(Math(txt_y))\n",
    "display(Math(txt_mu))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcae0fa8",
   "metadata": {},
   "source": [
    "Now that we have made the functional forms for each of our parameters, we'll sample all at once! As before, we need to define a overall log-likelihood function to pass into `emcee`. A large benefit of using the log-likelihood is that one can combine log-likelihoods linearly!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ddfcaa0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_Ltot(theta, nu, S, Sres, Serr):\n",
    "    T, y, mu = theta\n",
    "    T = T*u.K \n",
    "    return log_LBB(T, nu, S, Serr) + log_Ly(y, nu, Sres, Serr) + log_Lmu(mu, nu, Sres, Serr)\n",
    "\n",
    "pos_tot = [T_ini.value, y_ini, mu_ini] + 1e-6*np.random.randn(32, 3)\n",
    "nwalkers_tot, ndim_tot = pos_tot.shape\n",
    "\n",
    "sampler_tot = mc.EnsembleSampler(nwalkers_tot, ndim_tot, log_Ltot, args = (firas_freq, firas[\"monopole\"], firas[\"res\"], firas[\"sigma\"]))\n",
    "sampler_tot.run_mcmc(pos_tot, chain_len, progress = True);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb241064-bb28-46ea-9157-3ab02963d4ae",
   "metadata": {},
   "source": [
    "As before, let's plot the chains for each of the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d277bd71",
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = sampler_tot.get_chain()\n",
    "\n",
    "fig, axes = plt.subplots(sharex=True, nrows=3)\n",
    "labels=[r'$T_0$', r'$y$', r'$\\mu$']\n",
    "for i in range(3):\n",
    "    axes[i].plot(samples[50:, :, i], \"k\", alpha = 0.3)\n",
    "    axes[i].set_ylabel(labels[i])\n",
    "    axes[i].yaxis.set_label_coords(-0.1, 0.5)\n",
    "axes[2].set_xlabel(\"sample number\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cf5547a",
   "metadata": {},
   "source": [
    "Again, use the autocorrelation length to estimate the burnin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f18678a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = sampler_tot.get_chain()\n",
    "\n",
    "fig, axes = plt.subplots(sharex=True, nrows=3)\n",
    "labels=[r'$T_0$', r'$y$', r'$\\mu$']\n",
    "for i in range(3):\n",
    "    axes[i].plot()\n",
    "    axes[i].set_ylabel(labels[i])\n",
    "    axes[i].yaxis.set_label_coords(-0.1, 0.5)\n",
    "axes[2].set_xlabel(\"sample number\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36133312",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "flat_samples = sampler_tot.get_chain(discard = , thin = int(tau_y), flat = True)\n",
    "txt = \"\\mathrm{{{3}}} = {0:.7f}_{{-{1:.6f}}}^{{{2:.6f}}}\"\n",
    "labels = [\"T_0\", \"y\", \"\\mu\"]\n",
    "for i in range(ndim_tot):\n",
    "    mcmc = np.percentile(flat_samples[:, i], [16, 50, 84])\n",
    "    q = np.diff(mcmc)\n",
    "    txt_0 = txt.format(mcmc[1], q[0], q[1], labels[i])\n",
    "    display(Math(txt_0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f4d4b87-6fca-4894-84b7-a56b1d56661d",
   "metadata": {},
   "source": [
    "It is definitely useful to get an idea of how good of a fit our model is to our data, so we will define a simple sum of least-square function to evaluate. Luckily, the log-likelihood function that we've defined for the model evaluation is already written, and can be re-used to define our goodness-of-fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96984e4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "chisq = np.zeros((flat_samples.shape[0], 1))\n",
    "for i in range(flat_samples.shape[0]):\n",
    "    theta = [flat_samples[i, 0], flat_samples[i, 1], flat_samples[i, 2]]\n",
    "    chisq[i] = log_Ltot(theta, firas_freq, firas[\"monopole\"], firas[\"res\"], firas[\"sigma\"])*-2.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8843ad96-98fe-4550-9d23-91329b50265a",
   "metadata": {},
   "source": [
    "A good way to check how correlated or degenerate parameters are can be to plot corner plots ([python package](https://corner.readthedocs.io/en/latest/pages/quickstart/)). So let's give that a try here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b5d2a13",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(chisq.shape, flat_samples.shape)\n",
    "corner_data = np.concatenate((flat_samples[250:,:], chisq[250:,:]), axis = 1)\n",
    "corner(corner_data, show_titles = True, labels = [r\"T$_0$\", r\"$y$\", r\"$\\mu$\", r\"$\\chi^2$\"], title_fmt = \".3e\")\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4b11c23",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_mcmc = np.percentile(flat_samples[:, 1], 50)\n",
    "mu_mcmc = np.percentile(flat_samples[:, 2], 50)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(firas_freq, firas[\"res\"], \"k\", lw = 1)\n",
    "plt.plot(firas_freq, (func_mu(firas_freq, mu_mcmc)), color = 'xkcd:turquoise', ls = \"--\", label = \"fit (mu)\")\n",
    "plt.plot(firas_freq, (func_y(firas_freq, y_mcmc)), color = 'xkcd:orange', ls = \"-.\", label = \"fit (y)\")\n",
    "plt.ylim([-50, 50])\n",
    "plt.xlabel(r\"Frequency (s$^{-1}$)\")\n",
    "plt.ylabel(\"Residual intensity (kJy/sr)\")\n",
    "plt.title(\"Fits to residuals\")\n",
    "plt.tight_layout()\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc4da0dc",
   "metadata": {},
   "source": [
    "### Final consistency check\n",
    "\n",
    "We have not yet currently derived $y$-distortions or $\\mu$-distortions, so upper limits may be more appropriate to fit for. How would you report the results as upper limits rather than measurements centered around zeros?\n",
    "\n",
    "\n",
    "As a second point, the limits we derived are much tighter than the official limits. What are some potential sources of this difference? (There are many correct answers to this!)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3ac71f1",
   "metadata": {},
   "source": [
    "*I hypothesisize that...*"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}