Twinkle Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TWINKLE : THE SED GENERATOR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This tool, although a little rusty, will generate a stellar spectral energy distribution (SED) using empirical broad-band photometric data, fit to stellar models.\n",
    "\n",
    "#### Right now, the two free parameters are the stellar radius and temperature. The gravity and metallicty are held fixed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Input File Structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 JSON Parameter File"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The program works by feeding it a JSON parameter file that contains the names of folders and files that will be used by the code, boolean values that will dictate some fit hacks, plot commands and which photometry from the input file that will be used. \n",
    "\n",
    "For instance, the example JSON file is located and titled $\\textbf{./sed_paramfile.json}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1.1 File set up:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the JSON file, the \"folders\" section contains the following keys:\n",
    "\n",
    "    > \"topdir\" - Top directory name.\n",
    "    > \"codedir\" - directory under \"topdir\" that contains the code for twinkle.\n",
    "    > \"supportdir\" - Directory where the empirical standard colors, stellar models, and response functions are kept.\n",
    "\n",
    "and the \"files\" section contains the following keys:\n",
    "\n",
    "    > \"starfile\" - The file that contains empirical stellar data (photometry, parallax, etc).\n",
    "    > \"bv_colorfile\" - textfile that contains the empirical standard colors to use a reference.\n",
    "    \n",
    "If all the files and folders are located in a directory called $\\textbf{stinfo_topdir}$, then the following files should be structured in the following way:\n",
    "    \n",
    "    > topdir/codedir/starfile\n",
    "    > topdir/codedir/sed_paramfile.json\n",
    "    > topdir/codedir/plot_sed.py\n",
    "    > topdir/codedir/sed.py\n",
    "    \n",
    "    > topdir/codedir/supportdir/EMamajek_MSColors.txt\n",
    "    > topdir/codedir/supportdir/Models/\n",
    "    > topdir/codedir/supportdir/RSR/\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1.2 Photometric Bands:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use the JSON file to indicate which photometric band data to use for different purposes. There are  six main keys which are lists. Each key is an array of strings, where each string entry is the identifier for a specific photometric band.\n",
    "\n",
    "    > \"mags2use0_original\"\n",
    "        - All the bands to include at ANY POINT in the calculation.\n",
    "\n",
    "The rest are subsets of \"mags2use0_original\"\n",
    "\n",
    "    > \"mags4Phot0_original\"\n",
    "        - Bands that will be fit to the stellar model\n",
    "\n",
    "    > \"mags4scale0_original\"\n",
    "        - Bands that are used to calculate a first order constant to scale the raw model flux to observed         surface flux.\n",
    "\n",
    "    > \"scaleSEDbands\"\n",
    "        - Bands that will be used to re-scale the fitted photosphere. To be used with the boolean key \"scalephot\". \n",
    "    \n",
    "    > \"Remove_RedStars\"\n",
    "\n",
    "    > \"mags4Dust0\"\n",
    "        - Bands at which to calculate the excess flux and eventually one or two blackbodies.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example:\n",
    "    > \"mags2use0_original\" : [\"BJ\",\"VJ\",\"J2M\",\"H2M\",\"Ks2M\",\"W1\",\"W2\",\"W3\",\"W4\",\"HPACS160_flux\",\"HPACS100_flux\"],\n",
    "    \n",
    "    Will incorporate the Johnson B,V, 2MASS J,H,Ks, and WISE W1,W2,W3,W4, and Herschel 100, 160  photometric bands. The suffix \"_flux\" indicates that the input data is flux in jansky and not in units of magnitude.\n",
    "    \n",
    "    \n",
    "    > \"mags4Phot0_original\"  : [\"BJ\",\"VJ\",\"J2M\",\"H2M\",\"Ks2M\",\"W1\",\"W2\"],\n",
    "    \n",
    "    This indicates that only the Johnson B,V, 2MASS J,H,Ks, and WISE W1, W2 bands to use to fit the stellar photosphere. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### For all available filters, check out the $\\textbf{/Docs/available_filters.txt}$ file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 Stellar Photometric Data File"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The text file listed as \"topdir/codedir/starfile\" contains information for all the stars for which you want to create an SED. Each row contains information for a single star, where a column is the meta data for that star.\n",
    "\n",
    "The folder can be tab or space separated. \n",
    "\n",
    "The first row is the header. There are required and optional header names. They can be found in the file $\\textbf{/Docs/Starfile_ColumnDescription.txt}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can add as many unique photometric columns as you wish. As described in the $\\textbf{Starfile_ColumnDescription.txt}$, any column header key with the suffix $\\_flux$, should be in units of Jy. Otherwise, it will be assumed the units are in magnitudes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Running the Code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Set up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import sys, os\n",
    "import json\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from twinkle import twinkle\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### INITIATE INTERACTIVE MATPLOTLIB PLOTTING"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### If this doesn't work, change to \"%matplotlib inline\" and update both matplotlib and jupyter notebook to make nbagg backend to work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib nbagg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 Load JSON file\n",
    "##### Change path to where you store your paramfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "jfile = os.path.join(os.path.expanduser('~'),'Dropbox','Research','twinkle','twinkle','sed_paramfile.json')\n",
    "#jfile = 'Users/rpatel/sed_paramfile.json'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "script = open(jfile).read()\n",
    "specs = json.loads(script)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3  Select a star and feed JSON object into twinkle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "> /Users/rpatel/Dropbox (Personal)/Research/twinkle/twinkle/sed.py(1322)__init__()\n",
      "-> workingdir = directories.WorkingDir(specs['folders']['topdir'])\n",
      "(Pdb) c\n",
      "-------------------------------\n",
      "      Loading All Gridmodels   \n",
      "Loaded NextGen of g=40, met=0\n",
      "       Done Loading Models     \n",
      "-------------------------------\n",
      "Bands used to fit photosphere: [u'BJ' u'H2M' u'J2M' u'Ks2M' u'VJ' u'W1' u'W2']\n",
      "Bands used to scale photosphere: [u'BJ' u'H2M' u'J2M' u'Ks2M' u'VJ']\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "twinkle/utils/mosaic_tools.py:4131: VisibleDeprecationWarning: `rank` is deprecated; use the `ndim` attribute or function instead. To find the rank of a matrix see `numpy.linalg.matrix_rank`.\n",
      "  if np.rank(rr) != 2:\n",
      "twinkle/sed.py:765: RuntimeWarning: divide by zero encountered in log10\n",
      "  lam_arr_all, flux_arr_all = np.log10(lam_arr_all), np.log10(flux_arr_all)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi2 = 5.78\n",
      "Fitted Stellar Radius: 2.253 +/- 0.014 Rsun\n",
      "Fitted Stellar Temperature: 8923 +/- 46 K\n"
     ]
    }
   ],
   "source": [
    "# select a name in the stellar text file\n",
    "# check out the example file \"sample_stardata.txt\" as an example\n",
    "sname = 'Star_ID1'\n",
    "tauA = twinkle.Star(jfile, starname=sname)\n",
    "stoo = twinkle.STools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### AND THAT'S IT FOLKS! Your SED is lit!\n",
    "\n",
    "In the output you'll see which bands were used for the photospheric model fits, the associated chi2 value, and the optimal values for the stellar radius and temperature from the fit\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now Let's do some more, shall we?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. What Can I Access?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once a twinkle object has been instantiated, you can access various information for that star, and the functions in both twinkle.py and sed.py."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Star: Star_ID1\n",
      "distance: 64.1025641026\n",
      "NextGen 8923.73162413\n"
     ]
    }
   ],
   "source": [
    "print('Star:',tauA.starname)\n",
    "print('distance:',tauA.disti)\n",
    "\n",
    "# model used and fitted temperature (Kelvin)\n",
    "print(tauA.modeli, tauA.StarTemp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[u'BJ', u'H2M', u'HPACS100', u'HPACS160', u'J2M', u'Ks2M', u'VJ', u'W1', u'W2', u'W3', u'W4']\n"
     ]
    }
   ],
   "source": [
    "# All the bands included in the simulation.\n",
    "# fluxes are also calculated for each of these\n",
    "print(tauA.mags2use)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Plot Stellar SED"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### You can either plot the data with the built in \"plot_photosphere\" and \"plot_observedData\" functions, or access the global data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "/* Put everything inside the global mpl namespace */\n",
       "window.mpl = {};\n",
       "\n",
       "\n",
       "mpl.get_websocket_type = function() {\n",
       "    if (typeof(WebSocket) !== 'undefined') {\n",
       "        return WebSocket;\n",
       "    } else if (typeof(MozWebSocket) !== 'undefined') {\n",
       "        return MozWebSocket;\n",
       "    } else {\n",
       "        alert('Your browser does not have WebSocket support.' +\n",
       "              'Please try Chrome, Safari or Firefox ≥ 6. ' +\n",
       "              'Firefox 4 and 5 are also supported but you ' +\n",
       "              'have to enable WebSockets in about:config.');\n",
       "    };\n",
       "}\n",
       "\n",
       "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n",
       "    this.id = figure_id;\n",
       "\n",
       "    this.ws = websocket;\n",
       "\n",
       "    this.supports_binary = (this.ws.binaryType != undefined);\n",
       "\n",
       "    if (!this.supports_binary) {\n",
       "        var warnings = document.getElementById(\"mpl-warnings\");\n",
       "        if (warnings) {\n",
       "            warnings.style.display = 'block';\n",
       "            warnings.textContent = (\n",
       "                \"This browser does not support binary websocket messages. \" +\n",
       "                    \"Performance may be slow.\");\n",
       "        }\n",
       "    }\n",
       "\n",
       "    this.imageObj = new Image();\n",
       "\n",
       "    this.context = undefined;\n",
       "    this.message = undefined;\n",
       "    this.canvas = undefined;\n",
       "    this.rubberband_canvas = undefined;\n",
       "    this.rubberband_context = undefined;\n",
       "    this.format_dropdown = undefined;\n",
       "\n",
       "    this.image_mode = 'full';\n",
       "\n",
       "    this.root = $('<div/>');\n",
       "    this._root_extra_style(this.root)\n",
       "    this.root.attr('style', 'display: inline-block');\n",
       "\n",
       "    $(parent_element).append(this.root);\n",
       "\n",
       "    this._init_header(this);\n",
       "    this._init_canvas(this);\n",
       "    this._init_toolbar(this);\n",
       "\n",
       "    var fig = this;\n",
       "\n",
       "    this.waiting = false;\n",
       "\n",
       "    this.ws.onopen =  function () {\n",
       "            fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n",
       "            fig.send_message(\"send_image_mode\", {});\n",
       "            if (mpl.ratio != 1) {\n",
       "                fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n",
       "            }\n",
       "            fig.send_message(\"refresh\", {});\n",
       "        }\n",
       "\n",
       "    this.imageObj.onload = function() {\n",
       "            if (fig.image_mode == 'full') {\n",
       "                // Full images could contain transparency (where diff images\n",
       "                // almost always do), so we need to clear the canvas so that\n",
       "                // there is no ghosting.\n",
       "                fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n",
       "            }\n",
       "            fig.context.drawImage(fig.imageObj, 0, 0);\n",
       "        };\n",
       "\n",
       "    this.imageObj.onunload = function() {\n",
       "        fig.ws.close();\n",
       "    }\n",
       "\n",
       "    this.ws.onmessage = this._make_on_message_function(this);\n",
       "\n",
       "    this.ondownload = ondownload;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_header = function() {\n",
       "    var titlebar = $(\n",
       "        '<div class=\"ui-dialog-titlebar ui-widget-header ui-corner-all ' +\n",
       "        'ui-helper-clearfix\"/>');\n",
       "    var titletext = $(\n",
       "        '<div class=\"ui-dialog-title\" style=\"width: 100%; ' +\n",
       "        'text-align: center; padding: 3px;\"/>');\n",
       "    titlebar.append(titletext)\n",
       "    this.root.append(titlebar);\n",
       "    this.header = titletext[0];\n",
       "}\n",
       "\n",
       "\n",
       "\n",
       "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n",
       "\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n",
       "\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_canvas = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var canvas_div = $('<div/>');\n",
       "\n",
       "    canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n",
       "\n",
       "    function canvas_keyboard_event(event) {\n",
       "        return fig.key_event(event, event['data']);\n",
       "    }\n",
       "\n",
       "    canvas_div.keydown('key_press', canvas_keyboard_event);\n",
       "    canvas_div.keyup('key_release', canvas_keyboard_event);\n",
       "    this.canvas_div = canvas_div\n",
       "    this._canvas_extra_style(canvas_div)\n",
       "    this.root.append(canvas_div);\n",
       "\n",
       "    var canvas = $('<canvas/>');\n",
       "    canvas.addClass('mpl-canvas');\n",
       "    canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n",
       "\n",
       "    this.canvas = canvas[0];\n",
       "    this.context = canvas[0].getContext(\"2d\");\n",
       "\n",
       "    var backingStore = this.context.backingStorePixelRatio ||\n",
       "\tthis.context.webkitBackingStorePixelRatio ||\n",
       "\tthis.context.mozBackingStorePixelRatio ||\n",
       "\tthis.context.msBackingStorePixelRatio ||\n",
       "\tthis.context.oBackingStorePixelRatio ||\n",
       "\tthis.context.backingStorePixelRatio || 1;\n",
       "\n",
       "    mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n",
       "\n",
       "    var rubberband = $('<canvas/>');\n",
       "    rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n",
       "\n",
       "    var pass_mouse_events = true;\n",
       "\n",
       "    canvas_div.resizable({\n",
       "        start: function(event, ui) {\n",
       "            pass_mouse_events = false;\n",
       "        },\n",
       "        resize: function(event, ui) {\n",
       "            fig.request_resize(ui.size.width, ui.size.height);\n",
       "        },\n",
       "        stop: function(event, ui) {\n",
       "            pass_mouse_events = true;\n",
       "            fig.request_resize(ui.size.width, ui.size.height);\n",
       "        },\n",
       "    });\n",
       "\n",
       "    function mouse_event_fn(event) {\n",
       "        if (pass_mouse_events)\n",
       "            return fig.mouse_event(event, event['data']);\n",
       "    }\n",
       "\n",
       "    rubberband.mousedown('button_press', mouse_event_fn);\n",
       "    rubberband.mouseup('button_release', mouse_event_fn);\n",
       "    // Throttle sequential mouse events to 1 every 20ms.\n",
       "    rubberband.mousemove('motion_notify', mouse_event_fn);\n",
       "\n",
       "    rubberband.mouseenter('figure_enter', mouse_event_fn);\n",
       "    rubberband.mouseleave('figure_leave', mouse_event_fn);\n",
       "\n",
       "    canvas_div.on(\"wheel\", function (event) {\n",
       "        event = event.originalEvent;\n",
       "        event['data'] = 'scroll'\n",
       "        if (event.deltaY < 0) {\n",
       "            event.step = 1;\n",
       "        } else {\n",
       "            event.step = -1;\n",
       "        }\n",
       "        mouse_event_fn(event);\n",
       "    });\n",
       "\n",
       "    canvas_div.append(canvas);\n",
       "    canvas_div.append(rubberband);\n",
       "\n",
       "    this.rubberband = rubberband;\n",
       "    this.rubberband_canvas = rubberband[0];\n",
       "    this.rubberband_context = rubberband[0].getContext(\"2d\");\n",
       "    this.rubberband_context.strokeStyle = \"#000000\";\n",
       "\n",
       "    this._resize_canvas = function(width, height) {\n",
       "        // Keep the size of the canvas, canvas container, and rubber band\n",
       "        // canvas in synch.\n",
       "        canvas_div.css('width', width)\n",
       "        canvas_div.css('height', height)\n",
       "\n",
       "        canvas.attr('width', width * mpl.ratio);\n",
       "        canvas.attr('height', height * mpl.ratio);\n",
       "        canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n",
       "\n",
       "        rubberband.attr('width', width);\n",
       "        rubberband.attr('height', height);\n",
       "    }\n",
       "\n",
       "    // Set the figure to an initial 600x600px, this will subsequently be updated\n",
       "    // upon first draw.\n",
       "    this._resize_canvas(600, 600);\n",
       "\n",
       "    // Disable right mouse context menu.\n",
       "    $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n",
       "        return false;\n",
       "    });\n",
       "\n",
       "    function set_focus () {\n",
       "        canvas.focus();\n",
       "        canvas_div.focus();\n",
       "    }\n",
       "\n",
       "    window.setTimeout(set_focus, 100);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_toolbar = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var nav_element = $('<div/>')\n",
       "    nav_element.attr('style', 'width: 100%');\n",
       "    this.root.append(nav_element);\n",
       "\n",
       "    // Define a callback function for later on.\n",
       "    function toolbar_event(event) {\n",
       "        return fig.toolbar_button_onclick(event['data']);\n",
       "    }\n",
       "    function toolbar_mouse_event(event) {\n",
       "        return fig.toolbar_button_onmouseover(event['data']);\n",
       "    }\n",
       "\n",
       "    for(var toolbar_ind in mpl.toolbar_items) {\n",
       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
       "\n",
       "        if (!name) {\n",
       "            // put a spacer in here.\n",
       "            continue;\n",
       "        }\n",
       "        var button = $('<button/>');\n",
       "        button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' +\n",
       "                        'ui-button-icon-only');\n",
       "        button.attr('role', 'button');\n",
       "        button.attr('aria-disabled', 'false');\n",
       "        button.click(method_name, toolbar_event);\n",
       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
       "\n",
       "        var icon_img = $('<span/>');\n",
       "        icon_img.addClass('ui-button-icon-primary ui-icon');\n",
       "        icon_img.addClass(image);\n",
       "        icon_img.addClass('ui-corner-all');\n",
       "\n",
       "        var tooltip_span = $('<span/>');\n",
       "        tooltip_span.addClass('ui-button-text');\n",
       "        tooltip_span.html(tooltip);\n",
       "\n",
       "        button.append(icon_img);\n",
       "        button.append(tooltip_span);\n",
       "\n",
       "        nav_element.append(button);\n",
       "    }\n",
       "\n",
       "    var fmt_picker_span = $('<span/>');\n",
       "\n",
       "    var fmt_picker = $('<select/>');\n",
       "    fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content');\n",
       "    fmt_picker_span.append(fmt_picker);\n",
       "    nav_element.append(fmt_picker_span);\n",
       "    this.format_dropdown = fmt_picker[0];\n",
       "\n",
       "    for (var ind in mpl.extensions) {\n",
       "        var fmt = mpl.extensions[ind];\n",
       "        var option = $(\n",
       "            '<option/>', {selected: fmt === mpl.default_extension}).html(fmt);\n",
       "        fmt_picker.append(option)\n",
       "    }\n",
       "\n",
       "    // Add hover states to the ui-buttons\n",
       "    $( \".ui-button\" ).hover(\n",
       "        function() { $(this).addClass(\"ui-state-hover\");},\n",
       "        function() { $(this).removeClass(\"ui-state-hover\");}\n",
       "    );\n",
       "\n",
       "    var status_bar = $('<span class=\"mpl-message\"/>');\n",
       "    nav_element.append(status_bar);\n",
       "    this.message = status_bar[0];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.request_resize = function(x_pixels, y_pixels) {\n",
       "    // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n",
       "    // which will in turn request a refresh of the image.\n",
       "    this.send_message('resize', {'width': x_pixels, 'height': y_pixels});\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.send_message = function(type, properties) {\n",
       "    properties['type'] = type;\n",
       "    properties['figure_id'] = this.id;\n",
       "    this.ws.send(JSON.stringify(properties));\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.send_draw_message = function() {\n",
       "    if (!this.waiting) {\n",
       "        this.waiting = true;\n",
       "        this.ws.send(JSON.stringify({type: \"draw\", figure_id: this.id}));\n",
       "    }\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
       "    var format_dropdown = fig.format_dropdown;\n",
       "    var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n",
       "    fig.ondownload(fig, format);\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype.handle_resize = function(fig, msg) {\n",
       "    var size = msg['size'];\n",
       "    if (size[0] != fig.canvas.width || size[1] != fig.canvas.height) {\n",
       "        fig._resize_canvas(size[0], size[1]);\n",
       "        fig.send_message(\"refresh\", {});\n",
       "    };\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_rubberband = function(fig, msg) {\n",
       "    var x0 = msg['x0'] / mpl.ratio;\n",
       "    var y0 = (fig.canvas.height - msg['y0']) / mpl.ratio;\n",
       "    var x1 = msg['x1'] / mpl.ratio;\n",
       "    var y1 = (fig.canvas.height - msg['y1']) / mpl.ratio;\n",
       "    x0 = Math.floor(x0) + 0.5;\n",
       "    y0 = Math.floor(y0) + 0.5;\n",
       "    x1 = Math.floor(x1) + 0.5;\n",
       "    y1 = Math.floor(y1) + 0.5;\n",
       "    var min_x = Math.min(x0, x1);\n",
       "    var min_y = Math.min(y0, y1);\n",
       "    var width = Math.abs(x1 - x0);\n",
       "    var height = Math.abs(y1 - y0);\n",
       "\n",
       "    fig.rubberband_context.clearRect(\n",
       "        0, 0, fig.canvas.width, fig.canvas.height);\n",
       "\n",
       "    fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_figure_label = function(fig, msg) {\n",
       "    // Updates the figure title.\n",
       "    fig.header.textContent = msg['label'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_cursor = function(fig, msg) {\n",
       "    var cursor = msg['cursor'];\n",
       "    switch(cursor)\n",
       "    {\n",
       "    case 0:\n",
       "        cursor = 'pointer';\n",
       "        break;\n",
       "    case 1:\n",
       "        cursor = 'default';\n",
       "        break;\n",
       "    case 2:\n",
       "        cursor = 'crosshair';\n",
       "        break;\n",
       "    case 3:\n",
       "        cursor = 'move';\n",
       "        break;\n",
       "    }\n",
       "    fig.rubberband_canvas.style.cursor = cursor;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_message = function(fig, msg) {\n",
       "    fig.message.textContent = msg['message'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_draw = function(fig, msg) {\n",
       "    // Request the server to send over a new figure.\n",
       "    fig.send_draw_message();\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_image_mode = function(fig, msg) {\n",
       "    fig.image_mode = msg['mode'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.updated_canvas_event = function() {\n",
       "    // Called whenever the canvas gets updated.\n",
       "    this.send_message(\"ack\", {});\n",
       "}\n",
       "\n",
       "// A function to construct a web socket function for onmessage handling.\n",
       "// Called in the figure constructor.\n",
       "mpl.figure.prototype._make_on_message_function = function(fig) {\n",
       "    return function socket_on_message(evt) {\n",
       "        if (evt.data instanceof Blob) {\n",
       "            /* FIXME: We get \"Resource interpreted as Image but\n",
       "             * transferred with MIME type text/plain:\" errors on\n",
       "             * Chrome.  But how to set the MIME type?  It doesn't seem\n",
       "             * to be part of the websocket stream */\n",
       "            evt.data.type = \"image/png\";\n",
       "\n",
       "            /* Free the memory for the previous frames */\n",
       "            if (fig.imageObj.src) {\n",
       "                (window.URL || window.webkitURL).revokeObjectURL(\n",
       "                    fig.imageObj.src);\n",
       "            }\n",
       "\n",
       "            fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n",
       "                evt.data);\n",
       "            fig.updated_canvas_event();\n",
       "            fig.waiting = false;\n",
       "            return;\n",
       "        }\n",
       "        else if (typeof evt.data === 'string' && evt.data.slice(0, 21) == \"data:image/png;base64\") {\n",
       "            fig.imageObj.src = evt.data;\n",
       "            fig.updated_canvas_event();\n",
       "            fig.waiting = false;\n",
       "            return;\n",
       "        }\n",
       "\n",
       "        var msg = JSON.parse(evt.data);\n",
       "        var msg_type = msg['type'];\n",
       "\n",
       "        // Call the  \"handle_{type}\" callback, which takes\n",
       "        // the figure and JSON message as its only arguments.\n",
       "        try {\n",
       "            var callback = fig[\"handle_\" + msg_type];\n",
       "        } catch (e) {\n",
       "            console.log(\"No handler for the '\" + msg_type + \"' message type: \", msg);\n",
       "            return;\n",
       "        }\n",
       "\n",
       "        if (callback) {\n",
       "            try {\n",
       "                // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n",
       "                callback(fig, msg);\n",
       "            } catch (e) {\n",
       "                console.log(\"Exception inside the 'handler_\" + msg_type + \"' callback:\", e, e.stack, msg);\n",
       "            }\n",
       "        }\n",
       "    };\n",
       "}\n",
       "\n",
       "// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n",
       "mpl.findpos = function(e) {\n",
       "    //this section is from http://www.quirksmode.org/js/events_properties.html\n",
       "    var targ;\n",
       "    if (!e)\n",
       "        e = window.event;\n",
       "    if (e.target)\n",
       "        targ = e.target;\n",
       "    else if (e.srcElement)\n",
       "        targ = e.srcElement;\n",
       "    if (targ.nodeType == 3) // defeat Safari bug\n",
       "        targ = targ.parentNode;\n",
       "\n",
       "    // jQuery normalizes the pageX and pageY\n",
       "    // pageX,Y are the mouse positions relative to the document\n",
       "    // offset() returns the position of the element relative to the document\n",
       "    var x = e.pageX - $(targ).offset().left;\n",
       "    var y = e.pageY - $(targ).offset().top;\n",
       "\n",
       "    return {\"x\": x, \"y\": y};\n",
       "};\n",
       "\n",
       "/*\n",
       " * return a copy of an object with only non-object keys\n",
       " * we need this to avoid circular references\n",
       " * http://stackoverflow.com/a/24161582/3208463\n",
       " */\n",
       "function simpleKeys (original) {\n",
       "  return Object.keys(original).reduce(function (obj, key) {\n",
       "    if (typeof original[key] !== 'object')\n",
       "        obj[key] = original[key]\n",
       "    return obj;\n",
       "  }, {});\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.mouse_event = function(event, name) {\n",
       "    var canvas_pos = mpl.findpos(event)\n",
       "\n",
       "    if (name === 'button_press')\n",
       "    {\n",
       "        this.canvas.focus();\n",
       "        this.canvas_div.focus();\n",
       "    }\n",
       "\n",
       "    var x = canvas_pos.x * mpl.ratio;\n",
       "    var y = canvas_pos.y * mpl.ratio;\n",
       "\n",
       "    this.send_message(name, {x: x, y: y, button: event.button,\n",
       "                             step: event.step,\n",
       "                             guiEvent: simpleKeys(event)});\n",
       "\n",
       "    /* This prevents the web browser from automatically changing to\n",
       "     * the text insertion cursor when the button is pressed.  We want\n",
       "     * to control all of the cursor setting manually through the\n",
       "     * 'cursor' event from matplotlib */\n",
       "    event.preventDefault();\n",
       "    return false;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
       "    // Handle any extra behaviour associated with a key event\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.key_event = function(event, name) {\n",
       "\n",
       "    // Prevent repeat events\n",
       "    if (name == 'key_press')\n",
       "    {\n",
       "        if (event.which === this._key)\n",
       "            return;\n",
       "        else\n",
       "            this._key = event.which;\n",
       "    }\n",
       "    if (name == 'key_release')\n",
       "        this._key = null;\n",
       "\n",
       "    var value = '';\n",
       "    if (event.ctrlKey && event.which != 17)\n",
       "        value += \"ctrl+\";\n",
       "    if (event.altKey && event.which != 18)\n",
       "        value += \"alt+\";\n",
       "    if (event.shiftKey && event.which != 16)\n",
       "        value += \"shift+\";\n",
       "\n",
       "    value += 'k';\n",
       "    value += event.which.toString();\n",
       "\n",
       "    this._key_event_extra(event, name);\n",
       "\n",
       "    this.send_message(name, {key: value,\n",
       "                             guiEvent: simpleKeys(event)});\n",
       "    return false;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.toolbar_button_onclick = function(name) {\n",
       "    if (name == 'download') {\n",
       "        this.handle_save(this, null);\n",
       "    } else {\n",
       "        this.send_message(\"toolbar_button\", {name: name});\n",
       "    }\n",
       "};\n",
       "\n",
       "mpl.figure.prototype.toolbar_button_onmouseover = function(tooltip) {\n",
       "    this.message.textContent = tooltip;\n",
       "};\n",
       "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to  previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Pan axes with left mouse, zoom with right\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n",
       "\n",
       "mpl.extensions = [\"eps\", \"jpeg\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\"];\n",
       "\n",
       "mpl.default_extension = \"png\";var comm_websocket_adapter = function(comm) {\n",
       "    // Create a \"websocket\"-like object which calls the given IPython comm\n",
       "    // object with the appropriate methods. Currently this is a non binary\n",
       "    // socket, so there is still some room for performance tuning.\n",
       "    var ws = {};\n",
       "\n",
       "    ws.close = function() {\n",
       "        comm.close()\n",
       "    };\n",
       "    ws.send = function(m) {\n",
       "        //console.log('sending', m);\n",
       "        comm.send(m);\n",
       "    };\n",
       "    // Register the callback with on_msg.\n",
       "    comm.on_msg(function(msg) {\n",
       "        //console.log('receiving', msg['content']['data'], msg);\n",
       "        // Pass the mpl event to the overriden (by mpl) onmessage function.\n",
       "        ws.onmessage(msg['content']['data'])\n",
       "    });\n",
       "    return ws;\n",
       "}\n",
       "\n",
       "mpl.mpl_figure_comm = function(comm, msg) {\n",
       "    // This is the function which gets called when the mpl process\n",
       "    // starts-up an IPython Comm through the \"matplotlib\" channel.\n",
       "\n",
       "    var id = msg.content.data.id;\n",
       "    // Get hold of the div created by the display call when the Comm\n",
       "    // socket was opened in Python.\n",
       "    var element = $(\"#\" + id);\n",
       "    var ws_proxy = comm_websocket_adapter(comm)\n",
       "\n",
       "    function ondownload(figure, format) {\n",
       "        window.open(figure.imageObj.src);\n",
       "    }\n",
       "\n",
       "    var fig = new mpl.figure(id, ws_proxy,\n",
       "                           ondownload,\n",
       "                           element.get(0));\n",
       "\n",
       "    // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n",
       "    // web socket which is closed, not our websocket->open comm proxy.\n",
       "    ws_proxy.onopen();\n",
       "\n",
       "    fig.parent_element = element.get(0);\n",
       "    fig.cell_info = mpl.find_output_cell(\"<div id='\" + id + \"'></div>\");\n",
       "    if (!fig.cell_info) {\n",
       "        console.error(\"Failed to find cell for figure\", id, fig);\n",
       "        return;\n",
       "    }\n",
       "\n",
       "    var output_index = fig.cell_info[2]\n",
       "    var cell = fig.cell_info[0];\n",
       "\n",
       "};\n",
       "\n",
       "mpl.figure.prototype.handle_close = function(fig, msg) {\n",
       "    var width = fig.canvas.width/mpl.ratio\n",
       "    fig.root.unbind('remove')\n",
       "\n",
       "    // Update the output cell to use the data from the current canvas.\n",
       "    fig.push_to_output();\n",
       "    var dataURL = fig.canvas.toDataURL();\n",
       "    // Re-enable the keyboard manager in IPython - without this line, in FF,\n",
       "    // the notebook keyboard shortcuts fail.\n",
       "    IPython.keyboard_manager.enable()\n",
       "    $(fig.parent_element).html('<img src=\"' + dataURL + '\" width=\"' + width + '\">');\n",
       "    fig.close_ws(fig, msg);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.close_ws = function(fig, msg){\n",
       "    fig.send_message('closing', msg);\n",
       "    // fig.ws.close()\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.push_to_output = function(remove_interactive) {\n",
       "    // Turn the data on the canvas into data in the output cell.\n",
       "    var width = this.canvas.width/mpl.ratio\n",
       "    var dataURL = this.canvas.toDataURL();\n",
       "    this.cell_info[1]['text/html'] = '<img src=\"' + dataURL + '\" width=\"' + width + '\">';\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.updated_canvas_event = function() {\n",
       "    // Tell IPython that the notebook contents must change.\n",
       "    IPython.notebook.set_dirty(true);\n",
       "    this.send_message(\"ack\", {});\n",
       "    var fig = this;\n",
       "    // Wait a second, then push the new image to the DOM so\n",
       "    // that it is saved nicely (might be nice to debounce this).\n",
       "    setTimeout(function () { fig.push_to_output() }, 1000);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_toolbar = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var nav_element = $('<div/>')\n",
       "    nav_element.attr('style', 'width: 100%');\n",
       "    this.root.append(nav_element);\n",
       "\n",
       "    // Define a callback function for later on.\n",
       "    function toolbar_event(event) {\n",
       "        return fig.toolbar_button_onclick(event['data']);\n",
       "    }\n",
       "    function toolbar_mouse_event(event) {\n",
       "        return fig.toolbar_button_onmouseover(event['data']);\n",
       "    }\n",
       "\n",
       "    for(var toolbar_ind in mpl.toolbar_items){\n",
       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
       "\n",
       "        if (!name) { continue; };\n",
       "\n",
       "        var button = $('<button class=\"btn btn-default\" href=\"#\" title=\"' + name + '\"><i class=\"fa ' + image + ' fa-lg\"></i></button>');\n",
       "        button.click(method_name, toolbar_event);\n",
       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
       "        nav_element.append(button);\n",
       "    }\n",
       "\n",
       "    // Add the status bar.\n",
       "    var status_bar = $('<span class=\"mpl-message\" style=\"text-align:right; float: right;\"/>');\n",
       "    nav_element.append(status_bar);\n",
       "    this.message = status_bar[0];\n",
       "\n",
       "    // Add the close button to the window.\n",
       "    var buttongrp = $('<div class=\"btn-group inline pull-right\"></div>');\n",
       "    var button = $('<button class=\"btn btn-mini btn-primary\" href=\"#\" title=\"Stop Interaction\"><i class=\"fa fa-power-off icon-remove icon-large\"></i></button>');\n",
       "    button.click(function (evt) { fig.handle_close(fig, {}); } );\n",
       "    button.mouseover('Stop Interaction', toolbar_mouse_event);\n",
       "    buttongrp.append(button);\n",
       "    var titlebar = this.root.find($('.ui-dialog-titlebar'));\n",
       "    titlebar.prepend(buttongrp);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._root_extra_style = function(el){\n",
       "    var fig = this\n",
       "    el.on(\"remove\", function(){\n",
       "\tfig.close_ws(fig, {});\n",
       "    });\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._canvas_extra_style = function(el){\n",
       "    // this is important to make the div 'focusable\n",
       "    el.attr('tabindex', 0)\n",
       "    // reach out to IPython and tell the keyboard manager to turn it's self\n",
       "    // off when our div gets focus\n",
       "\n",
       "    // location in version 3\n",
       "    if (IPython.notebook.keyboard_manager) {\n",
       "        IPython.notebook.keyboard_manager.register_events(el);\n",
       "    }\n",
       "    else {\n",
       "        // location in version 2\n",
       "        IPython.keyboard_manager.register_events(el);\n",
       "    }\n",
       "\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
       "    var manager = IPython.notebook.keyboard_manager;\n",
       "    if (!manager)\n",
       "        manager = IPython.keyboard_manager;\n",
       "\n",
       "    // Check for shift+enter\n",
       "    if (event.shiftKey && event.which == 13) {\n",
       "        this.canvas_div.blur();\n",
       "        event.shiftKey = false;\n",
       "        // Send a \"J\" for go to next cell\n",
       "        event.which = 74;\n",
       "        event.keyCode = 74;\n",
       "        manager.command_mode();\n",
       "        manager.handle_keydown(event);\n",
       "    }\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
       "    fig.ondownload(fig, null);\n",
       "}\n",
       "\n",
       "\n",
       "mpl.find_output_cell = function(html_output) {\n",
       "    // Return the cell and output element which can be found *uniquely* in the notebook.\n",
       "    // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n",
       "    // IPython event is triggered only after the cells have been serialised, which for\n",
       "    // our purposes (turning an active figure into a static one), is too late.\n",
       "    var cells = IPython.notebook.get_cells();\n",
       "    var ncells = cells.length;\n",
       "    for (var i=0; i<ncells; i++) {\n",
       "        var cell = cells[i];\n",
       "        if (cell.cell_type === 'code'){\n",
       "            for (var j=0; j<cell.output_area.outputs.length; j++) {\n",
       "                var data = cell.output_area.outputs[j];\n",
       "                if (data.data) {\n",
       "                    // IPython >= 3 moved mimebundle to data attribute of output\n",
       "                    data = data.data;\n",
       "                }\n",
       "                if (data['text/html'] == html_output) {\n",
       "                    return [cell, data, j];\n",
       "                }\n",
       "            }\n",
       "        }\n",
       "    }\n",
       "}\n",
       "\n",
       "// Register the function which deals with the matplotlib target/channel.\n",
       "// The kernel may be null if the page has been refreshed.\n",
       "if (IPython.notebook.kernel != null) {\n",
       "    IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n",
       "}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAYAAAB91L6VAAAgAElEQVR4XuydCZyN1f/H33ZZSyotKnubQpSlEClFpKS0KaWi5Zc2S8uvn5SIFkX1J0WkVBQVoqwpCZUlyVYq0mZfBuP/+t65M8YYc+/Mfe69z/I5r5cX7jzPOd/z/j4z7znPc55z8qEiAiIgAiIgAiKQcAL5Et6iGhQBERABERABEUAC1kUgAiIgAiIgAkkgIAEnAbqaFAEREAEREAEJWNeACIiACIiACCSBgAScBOhqUgREQAREQAQkYF0DIiACIiACIpAEAhJwEqCrSREQAREQARGQgHUNiIAIiIAIiEASCEjASYCuJkVABERABERAAtY1IAIiIAIiIAJJICABJwG6mhQBERABERABCVjXgAiIgAiIgAgkgYAEnAToalIEREAEREAEJGBdAyIgAiIgAiKQBAIScBKgq0kREAEREAERkIB1DYiACIiACIhAEghIwEmAriZFQAREQAREQALWNSACIiACIiACSSAgAScBupoUAREQAREQAQlY14AIiIAIiIAIJIGABJwE6GpSBERABERABCRgXQMiIAIiIAIikAQCEnASoKtJERABERABEZCAdQ2IgAiIgAiIQBIISMBJgK4mRUAEREAEREAC1jUgAiIgAiIgAkkgIAEnAbqaFAEREAEREAEJWNeACIiACIiACCSBgAScBOhqUgREQAREQAQkYF0DIiACIiACIpAEAhJwEqCrSREQAREQARGQgHUNiIAIiIAIiEASCEjASYCuJkVABERABERAAk7gNVC4cOF9Rx11VAJbVFMiIAIiIAJOEvjtt99SgCJO1CkBO0ExyjqOP/74fb/++muUR+swERABERABtxHIly/fb8AJTsQlATtBMco6JOAoQekwERABEXApAQnYpYmJFJYEHImQvi4CIiAC7iYgAbs7P4eMTgL2aOIUtgiIgAiECUjAHr0UJGCPJk5hi4AIiIAE7O1rQAL2dv4UvQiIgAhoBOzRa0AC9mjiFLYIiIAIaATs7WtAAvZ2/hS9CIiACGgE7NFrQAL2aOIUtgiIgAhoBOzta0AC9nb+FL0IiIAIaATs0WtAAvZo4hS2CIiACGgE7O1rQAL2dv4UvQiIgAhoBOzRa0AC9mjiFLYIiIAIaATs7WtAAvZ2/hS9CIiACGgE7NFrQAL2aOIUtgiIgAhoBOzta8DNAh4yZCK33XYpDRo0ZNCglzjrrOrehq3oRUAERCAOBDQCjgPURFTpVgHv2bOHQoUKHYBg+PBplC/fmMqVoXz5RNBRGyIgAiLgfgISsPtzlG2EbhXw22+/Tfv27TPFfAsUfA5OeZ/8y9vxykvF6NTJo9AVtgiIgAg4SEACdhBmIqtyq4B79XqR//73njCKE6DkLOhwMZRdDn9VheFT+GXJiRoJJ/JiUVsiIAKuJCABuzItkYNyq4Dz5cuXKfjGcMKjcGvT/Z8N/ZyLTltK794NqFOnRujztWthxQp0izpy2nWECIiAjwhIwB5NpjcEfAKwEh4vsp/y47uASsCvFCxYiMGDU+jSBQoXhpQUGDwY3aL26DWpsEVABHJHQALOHS/XHO0NARuuW6DEY5BSBArvgq29gNf236JmFRTcC6eMg2VtKEhRVq3SZC3XXGgKRAREIG4EJOC4oY1vxW4U8L59+8ifP382HbeRcCVKl97Gpk3fHHiLuuT/QYeWGc+Ii74zjYljjqNx4/jyU+0iIAIikGwCEnCyM5DH9t0o4L1791KwYMFD9uibb77h9ts7M3/+PP74YwPHHFMLTngLbm2YcU6B179g9az6mqSVx+tCp4mACHiHgATsnVwdEKkbBfzNNwupU6fWIYnaCLlFi1v55JPXsH+HJmzluw3++3/7Bdy7Mx9/dB0XX9zAo5lR2CIgAiIQHQEJODpO8TyqOdA93IAtGWVThr+N1KAbBTxkyOvcdlvHbEP//ff1HHvsMVx8cUc+/fT1DAHXq3cpXy5ayOjh82jfoQ7XX34NI0c+x6OPvkCvXumvM0Wioa+LgAiIgPcISMDuyZm9v/M1cA6wL1JYbhTwoEGvcNddnQ8IvUCBQuzduzskXCs33NCJkSOHZgj45Zdf45lnhrJy5Rzat7+X0aOfp3Hj65k7dxKXXNKSN98cRPHixSPh0NdFQAREwHMEJGD3pMzEezNwoMEOEZ8bBXzgO8BpgU+Y8BGXXdYyQ8BXXHEH48a9miHgDRs2cNRRRx3UywoVTmPNmh9Cn0+e/CUXXVTXPZlSJCIgAiLgAAEJOG8QBwKtgJOAmlluGVcBhgNlgU3ATcCSKJp5EpgFTIriWLwi4I8++piWLVtkCPjyyzvz4YevZAh406ZNlCpV6qAuf/zxLFq23D856+KL2zJp0rvRoNExIiACIuAJAhJw3tJkZlgFzAYuzyLgz4ERwBtAW6AbUAc4DRicpTmT7dPhz+YA9vJNSjQheUXAEydO5JJLLskQcIcOnRkxIk3ANmu6QIEC2XZ36tTpNGt2QehrZ59dn/nz59CgwaVMmfIehx12WDSIdIwIiIAIuJqABBxbetZkEfDRwAqgDLAHsOe664Dzwp8fqrWKYRG3yyGc+wD7EyqlS5c+fuPGjbFFn4uzd+7cSZ8+z9OrV49cnAWTJk2iefPmGQK+8srOjB2bJuCcyowZs2jcuCHNmt3OO+/0oUSJErRrdwsffPAmM2bMpWFDu2OvIgIiIALeJSABx5a7rAI+G3gLqJapWptYZbOcbWR8qNIV2ACMijYcp0fAGa8FAevX/8mGDetJTYWaNc+MNqSDjrP3fm1rwrPOOivbW9A5Vfzll19Sv379A0RtMV5ySUcmT36D//3v/3jsMW2rlOfk6EQREIGkE5CAY0uBUwLOdRSxCnjYsBHcckuHg9odNmwUHTtel+t4sjshu1HuDTfcwciRaZOwcivg9OPTJ3u1aHED11xzJaNGjefddwdSooRmSzuSOFUiAiKQEAIScGyYnboFnesoYhVwdjOWcx1EhBOyk2ybNrfzwQf/F1HAc+bYM98G2R53qNgnT57KRRdl2nnJ6Q6pPhEQARFwkIAEHBvMrAK22qaHJ2ClT8Ky28+1Y2vm4LO9KuDWrW9n/PjIAv7iiy8477zzshXw+vXrOeKIIyhatOhBYHr06MWTTz6StsqWigiIgAi4mIAEnLfkvAq0AMoBfwNbgMrhquz5r8n3SGBz+N3eRXlr5tBnuVnA/fo9y0MP3ZetPK+//nZGjYpNwOlUHnusN0888Wi2kCLd4nY6H6pPBERABHJLQALOLTGXHO9mAZv8Mk/qyoysdevbGD9+SMRb0LNnz+b888+PeFxOI93t23dw2GEHj5JdkkKFIQIiEHACErBHLwC3C/hQWFu1uo0JE5wT8NChr9Op0/71pwsXLkpKys5Q87YM5ooVP3HyybZeiooIiIAIuIuABOyufEQdjRsFXKtWIxYsmJHjqPXaa29j9OjIAp41axYNGzaMOAI2YMOGjeWWW64MsStUqDi7d2/j8cf78fjjD4U+27p1q9aTjvrK0oEiIAKJIiABJ4q0w+3EKuDKlU9n5cqluY7K3u2tXTv7OWXt2nVgzJjhOUoz2hHwzJkzadSoUVQCtk7UrNmab78dz7HHVmPduh95/PH+PP74A+GRcGG++upbatc+Ndf91QkiIAIiEC8CEnC8yMa53lgFfNddDzBo0IBcR7lu3TqOPfbYg85bvny5rU9NsWLFcqxzw4a/Wb78N847L+cFPnIrYGt006bN7Nmzl7Jly/DYYwPo1ev+UCzNm1/LpElvMXTo21x00dWsWAGVK0P58rnuvk4QAREQAccISMCOoUxsRbEK+IMPPqZNm5a5Dtp2Lzr6aFtxc39Zt2495codk+u6cjphwYIFnH322VGPgDPXZROz3nvvQ9q2bR36+Ouv59Gz5wCmTi1OvkIDKVj9Q1KXtOHlF4vSSYtpOZo3VSYCIhA9AQk4elauOjJWAa9evZqKFW0J6tyV7GY3x+uVn3///Tf0vm9eS/oM6Xnz5lG2bE0q1vidfTdcCGWXw19VKTByGqu/P04j4bwC1nkiIAIxEZCAY8KXvJODIOBY6aYLeO7cuWzbdg6X3PYlu66vn1Ft4TfnMHloPRrbHlQqIiACIpBgAhJwgoE71VwyBXzddT146630XRTJ021ipzjkVE+JEkeybds/2J7DmzaVokLFfex9JH/GKfl67eHnNQU0Ak5EMtSGCIjAQQQkYI9eFMkUsCHLvABGvG5Bx5oa2284f/78GbEOGQKdH1rH3h0FoNBOSvA6Tz99HHfeqQfBsbLW+SIgArknIAHnnpkrzohVwKtWraJSpUq57ku6bL0g4Ow6t3Yt1Kp1O3/9NZm//15AvXrNadmyAwMG3JlrFjpBBERABGIhIAHHQi+J57pFwD/+uJKqVXM/mSuJ6Khe/WyWLfue3bt3888/Gzn99Itp164ZL7zQO5lhqW0REIGAEZCAPZpwtwjYrbefc0rriSfW4rfflrB3767QYfYa1XHHHctll93N+PEDPXpFKGwREAGvEZCAvZaxcLzJEPBpp53JkiXfhSJIvwXtRQEff3xN/vjjB/bsSVsz2srLL4+lS5cruemm+3j99QHYrWot2OHRbw6FLQIeISABeyRRWcOMVcArV66ksi0HlYuSWbZFipQkJWWra2dA59Stu+9+kLVr1/LBB28fcNiFF17JZ5+NBW6hQJEXyX/6B1qwIxfXhw4VARHIHQEJOHe8XHN0sgW8aNEizjzzTE8KOKck5stXHkrOgg4Xa8EO11ztCkQE/ElAAvZoXp0WcIUKZ7F6ddrt5UOVzCPgJUuWcMYZZ/hOwH37fkP3lzbDrU0zMGjBDo9+kyhsEXA5AQnY5Qk6VHhOC/issy7ku++mHtTc3Lnfce65Z4U+zyzgX375hZNOOsl3ArZnvyeetBv+WziDRf4n9rJmdX4t2OHR7xWFLQJuJSABuzUzEeKKVcArVqygSpUqGa107foUzz3X86BWM6/9nHXC1fbt2yPufuRFvIMG7eKu7n9AShEovAu29mLbtoG+7KsX86OYRcAvBCRgj2bSaQFnt8lC+qjXyzOe85ren37ayauvTmXAgM7Ar6Fqdu7cSZEiRfJapc4TAREQgQMISMAevSAk4Pgnbs2aNVSrdhYpKZszGrOFO444onT8G1cLIiACvicgAXs0xRJw4hKXedlNa3XDhr846qgjExeAWhIBEfAlAQnYo2mNVcA//fQTVatWzei9bkEf+kLo1KkvQ4d2P+CAjRs3seGfwgz/ehwdzmlDlQpFPXolKWwREIFkEZCAk0U+xnYl4BgB5vL0AgUKkJqauv+sUtdDh7lw5E/wd1X6nTaNB+84Lpe16nAREIEgE5CAPZr9RAi4RImj2bLlD55/fjBdu97pu1eOcpN629bQ7hKcc851jB79OpUaz4Vbzs+oosDrc1g9q55eVcoNVB0rAgEnIAF79AJwUsCffvo1zZrVOWCPX8Ny1FEnsWHDGjZv3kzp0qUDLWB7DnzSSXVZs+ZLpk2DFi1T2fFQgYyrp2jfVCZ+ko/GjT16QSlsERCBhBOQgBOO3JkGYxXw8uXLqVatWiiY7Pb4tc+PPPJE/vrrZ7Zs2UKpUqUCLeAyZcpRvXodZsyYENqooWJF2FN0HaSUgMJbKbizHKtW5dMI2JnLW7WIQCAISMAeTXM8BHzbbXczZMhLGUSqVKnP8uVfSMDA33//TeHChSlZsmSIz5AhcNttuylUKB+7d0Pt2h8wb15bj15NClsERCAZBCTgZFB3oM14CLhr1548/3yfjOi2bNlKiRLF2bp1a0g8Xtx60AHUh6zCRsIrV9qdgk2ceebhvPrqm9x22/XxbFJ1i4AI+IiABOzRZMZDwHfd1ZNBg/YLOF2427Zto0SJEhJwDtfKunV/cOqpjejR4366devk0atKYYuACCSSgAScSNoOthWrgH/88UdOOeWUUETpor3//h48++zTGVGmf25rPhcvXlwCjpC/adO+oUmTOnTp0o1Bg/ZzdDDtqkoERMBHBCRgjyYzHgL+z396MHCgBBzLJdGhw4OMGNGfwYM/oHPn1rFUpXNFQAR8TkAC9miCYxVwdrOg77yzB4MHHyzgHTt2hHYC0jPg6C6WM89sxKJFM3nooWfp27drdCfpKBEQgcARkIA9mvJYBWwytcUlMt+C7tHjvzz9dK+DbkHbLkCHHXaYBBzltbJgwULOPrtW6OhRoz7i2mtbRHmmDhMBEQgSAQnYo9mOVcDW7azbDD72WH+eeOJBCTjGa+Lff/+lTJkylChRiq1bN7N48RJOP/20GGvV6SIgAn4jIAF7NKNOCbhs2VP588+lIQo9e/anT5+DBbxr1y6KFi2qEXAur5X77nuRTz4ZzY8/fslnn82nSZO0UbGKCIiACIQHQb8BJzhBI58TlaiO6Ag4IeAiRUpxxBHlWb9+SajRJ598hkceeeigEXBKSkpoI3o9A44uN5mP2rRpE4cffjilS1dlxYo5lC2rbQxzT1FniIA/CWgE7NG8OiFgu0V63HHlWb48TcCPPPIMTz4pATt9SXz11XdMnz6dHj3uZeXKX6lY8Xinm1B9IiACHiQgAXswaRayEwIuXLg0hx9+Ahs2pAm4e/d+9O3b7aAR8O7du0PLMGoEnPeL5fPPp9G0aRPKlDmDZcumsXFrCe0lnHecOlMEfEFAAvZoGp0QcNGipTniiBNYty79FnQ/HnkkTcBff/09depUD/17z549FCpUSAKO4VpZv349l132H9q2bcTQMf+wosEI7SUcA0+dKgJ+ICABezSLTgi4ZMnSlCt3Aj/9lCbgnj370adPmoAzj3b37t1LwYIFJWAHrpVfftnHSQ1mw60NM2rTXsIOgFUVIuBBAhKwB5Pm1C3ocuVqYNvsLV066aBb0BJwfC4M20u46YU72PdYsYwGtJdwfFirVhFwOwEJ2O0ZOkR8ToyAt27dFlqMo1ixw0Kt9O7dl0cf7X7QCDg1NZUCBQpoBOzAtZK2l/A+9hRdn2kv4WNZtQrtJewAX1UhAl4iIAF7KVuZYnVCwFm73r17X/r2PVjA6atmaRKWMxeL7SXcpQsUKrSPHTv20qbNV4wde54zlasWERABzxCQgD2TqgMDjYeAu3V7mn79ehw0ApaAnb9I0vcSLlXqH9q2bcB99/Xhrrsud74h1SgCIuBaAhKwa1OTc2DxEPATTzzNY48dLGCLxJat1Ag4PhfLmDGTuPrqS3juuaHce+8t8WlEtYqACLiOgFcE/LlD5PYBTR2qK6nVxEPA3br1pV+/g29BS8DxTXXmnamGDfuAm2/WNobxJa7aRcAdBLwi4FSHcJmACzhUV1KriY+As78FLQHHN9V//PEH5cqVy2hk/vzvqFXrzPg2qtpFQASSTsBLArZ3ZfrGQMyGdhdJwIcmqFvQMVxdDp2avkOVSfnoo492qFZVIwIi4EYCXhLwG0DHGCC+DtzoQgEXsm1jAftp+zHwTDR9jMcI+MEH+9C/f89Q81mf9+oZcDRZif2Yiy5qz5Qpb4cq+u23dRx33P6Rcey1qwYREAE3EfCKgJcB7wMPxwCvD2DTTE+NoY54nHo1cCzwPGA/ee8G/ozUUDwE/NBDfXjmmewF/OWX31KvXo1IYenrMRKwd7Nfeuld3nprBOvW/cG0WV/w9ncT6XBOG6pUKBpj7TpdBETATQS8ImA3MXM6Frs1Pg/4DPhv+N+fRGokHgKeOfMLGjVKex9VM54jZSD+X89Xqj10mA9lf4K/q9LvtGk8eMdx8W9YLYiACCSEgAScN8wDgVbASUBN4NtM1VQBhgNlgU3ATUDaYsvZlzbhOky+U4DXgNGRwoqHgK3NF14YSffuXdmxI+IgPFKI+noMBOw94ZPPm0NqxwYZtWjN6BiA6lQRcCEBCThvSbGV9FcBs8O3tTML2F6ZGgHYM+u2gO1uUAc4DRicpTmbWGbPfF8CTNxmvaHh0XCOkcVLwHnDobOcJmBrRrdouY8dD+XPqFprRjtNWfWJQHIJSMCx8V+TRcA2kWoFUMZ28bP1K4B1gN3Xtc9zKnasjZzvALZnc+B9gP0JldKlSx+/cePG2KLX2a4lkLZmNOwpui5jzej8249hzZr8WjPatVlTYCKQOwJuF/BRgD30WglszaZrJcO3b2fmrtuOHZ1VwGcDbwHVMrXwte11DxxqMZETwyNme9fZJmKNjyY6jYCjoeTtY/avGQ27d++jSJFePPnkGfznP1d6u2OKXgREIETArQIuCAwJvzZkcaYArwK2TuKOTLk7F5iTxFeLnBBwni5FCThP2Dx3Uvqa0ZUqQUrKrzRs2IBJkyZRvbrbJvN7Dq0CFoGkE3CrgO1Wa2/gSWA+0Ai4B1gONAf+CJNzm4BjuQWdq4tBAs4VLt8cXKlSXVatmsuXXy6mbt3TfdMvdUQEgkjArQJeHF6cwt7dTS82iemD8Gj34vAzVbcJ2GKdHp6AlT4Jy24/13b64pKAnSbqjfo++mgOr776FjNmTObrr6dwyikneyNwRSkCInAQAbcKeBtwKTAjS8RHAjZz2J6bXgLYKlLJuAVtt8NbALZM0d/AFqByOFZ7/mvytVg3AzcDi5y+9iRgp4l6q7433viI229vG1otq2zZI7wVvKIVAREIEXCrgH8GHgTGZJOnEuGJSjbhydaGfiKJz4CTdhlJwElD75qGL7vsPpYvX8rSpR9ToIAv9hhxDVsFIgKJIOBWAY8FdgLXHgJCYeAdwPZt880OR7lJuAScG1r+PNZWKzv66BOoUqUGX3zxUWjPZhUREAHvEHCrgG0Bi/uBluFbvNkRtRUKbGELex5cwTvInYlUAnaGo9drMQnnz58/JN+9e/dKwl5PqOIPFAG3CjhQSchLZyXgvFDz5zknnngWa9d+z3/+8wTPP/+IPzupXomADwlIwB5NqgTs0cTFKWzbRalixWrccks3+vSxDbVUREAE3E7AawK+ApgM2CzpQBcJONDpz7bztWq1YeHCD3jrram0b99UgERABFxOwGsC3hvez9cW5Ah0kYADnf5sO79nzx4WLlxOkyaXMWXKeC3UoUtEBFxOwGsCtvWSTwmviOVytPENTwKOL18v1/7uu9Np1+4C1qz5nZNOOtbLXVHsIuBrAhKwR9MrAXs0cQkK22ZF2+tJCxbMokQJe3VeRQREwG0EJGC3ZSTKeCTgKEEF+LCePZ9h9OiJLFs2kSJFigSYhLouAu4kIAG7My8Ro5KAIyLSAUCbNt3Zvn0tkyaN1DvCuiJEwGUEJGCXJSTacCTgaEkF+7jU1NTQMpWnn16TxYsXBBuGei8CLiMgAbssIdGGIwFHS0rHpS9R2a3bKzz99O0CIgIi4BICErBLEpHbMCTg3BIL7vEtWtxD586tuPXW+xk8+HGuuKJNcGGo5yLgIgISsIuSkZtQJODc0NKxRmDt2rWceOKJjBnzKVdd1UxQREAEkkzAawLuANhOSbb/bqCLBBzo9Oe5808/PYY+fe5h4cI5VKxYMc/16EQREIHYCXhNwLH32Cc1SMA+SWQSuvHxxxO56ab7WLFiLqVLl0pCBGpSBETACEjAHr0OJGCPJs4lYadPzNq+fTuHHXaYS6JSGCIQLAJeF7CtLmBr7f0ZtA0aJOBgfaM63dsXXvg/7r33dho2vJLPP38n9KqSigiIQGIJeFXAhwNDgMuB/GFkq4H5wDfhv+3fmxKLM3GtScCJY+3Hlmzjhl9++Z127XpQo8aZDB3azY/dVJ9EwNUEvCrgoUBHYB+wHrAHWcXDpO2z9LIKmAdc6+os5CE4CTgP0HTKQQQ2btzMEUeU5p57nuOFF+4VIREQgQQS8KqAfwdKA43CI15DVhU4G6gV/rtm+BgTsu/ur0nACfwu8XlTK1asom7dNgwY0IcOHS71eW/VPRFwDwGvCngz8JktdRsBZaWwjMe4B7kzkUjAznBULWkEli5dxSWXXMygQSNo2bKesIiACCSAgFcF/GX41nMkAScAYXKakICTw93Prc6c+S0tW7Zl0aJZGfsI79yzk3HLxtHmlDYULVjUz91X30Qg4QS8KmBb0LYPcCKwNeHUXNCgBOyCJPgwhN69h9Cv30OsWbOaHQW3UXFgRVL2plC1TFWm3TSN40oe58Neq0sikBwCXhVwPmAGsCE8wSolOfiS16oEnDz2fm7Z3guuXPk8jjiiPIMn3EfjNxtndHdOxznUK6/b037Ov/qWWAJeFbBRsklYX4eXpXwCmAgERsQScGK/UYLU2r59+7j99nuZ+3UNvm9jLxuklVePTeW22+x3XxUREAEnCHhVwCeHhWszn+0ngs103gl8F37tyN4FttePljkByY11SMBuzIp/YlqzZi8VKuyDEn9CSgkovJWCO49l1SooX94//VRPRCCZBLwq4HFAa2AHsBgoCVTJ9LpR+rvA9nzYdiG/IJmQ49G2BBwPqqozncC0aXDppans3Jm+zg3YipWffAKN99+VFjAREIEYCHhVwP+ER721AVsBy0oxwN79tXeB7XP7kz5C1nvAMVwkOjV4BNauBdssac+e/X0vWBCNgIN3KajHcSTgVQFvBKYAV0VgUyIs5VlxZJiUqjUCTgr2QDU6ZAh06QIFC9pIeC9PPbWJHj3KBoqBOisC8STgVQFPA7YDLeIJx811S8Buzo5/YrOR8MqV8Nln43nxxQdZuHAaFSroVST/ZFg9SSYBrwq4HTAMsJWu/kgmwGS1LQEni3xw233iiRGMHDmAuXNncvjh9hKCigiIQCwEvCpg++5/EagAtAL+jQWCF8+VgL2YNe/H/PDDLzB69ASWLfuEwoULe79D6oEIJJGAVwW8F/gNOB74FXgWGAusTSLLhDYtAScUtxrLROCSS+6lWLFNvPfeMPLl03vBujhEIK8EvCpgm4RlWxCml/TXjmxrQnsHOP09YPv7r7zCcfN5ErCbs3D/io4AACAASURBVOPv2Pbu3UvTpm2pWPEchg3r4e/OqnciEEcCXhWwIUnf6Sh9+0H7+4gwq8x7Av8SvlUdR4yJr1oCTjxztbifwJYtW6lcuRG9ej3E7bdfLTQiIAJ5IOBlAWfXXVshy94DTt8X2N4FNinrPeA8XBw6RQRyIvD77+upW7chgwcP1xaGulREIA8E/Cbg7BDYjkk2CvZV0QjYV+n0bGdmzfqOSy+9gnnzPuOUU+z3XxUREIFoCQRBwNGy8NRxErCn0uXrYIcPn8DTTz/CV1/NonTpzFMzfN1tdU4EYibgVQHbFi03Af8DPouZggcrkIA9mDQfh9yly7N89tlUliwZT0Fbs1JFBEQgIgGvCti2HrQNFk7IYZaz/RTItJJtRBaeOkAC9lS6fB+sbWFYv/7lbNjwMytWLAzU60k79+xk3LJxtDmlDUULFvV9rtVB5wh4VcD2vu8qoFEOKGyd6P7A3cB455C5oyYJ2B15UBT7Caxdu5YzzjiX5s1v4Z13bItu/5ffNv9Gk+FNWP7PcqqWqcq0m6ZxXEkt1en/zDvTQ68K2Pb+fRe4IQcM9quoLVM5GbClK31VJGBfpdM3ndm+fQfnn9+M9u078MADnXzTr0N15Mu1X1J/WP2ML8/pOId65ev5vt/qoDMEvCrgDcBc4LIIGGYD5YDKzuByTy0SsHtyoUgOJPDnn39TtWpDBg0ayLXXNs0Wj19u29qt9/y99u+ZnPpYaqBuv+vaj42AVwVsWxHa3r92ryclBwSjgZZAydgwue9sCdh9OVFE+wksWPAjLVpcwuzZU6lUqeIBaPx223bdlnWUKFyCrSlbObbksboMRCBqAl4VcBfgJeAR4KkcejsVsPtDxaIm4pEDJWCPJCrAYb733sd07vwwy5ZN48gj0xepA922DfBFoa4fQMCrAj4M+D68xKRJ+Ols8noMsCK8QcNpfsu7BOy3jPqzP/fd9xKTJ7/DwoWfZeyepNu2/sy1epV7Al4VsPX0TMBGuEeGZfwK8AWwBTDh9gZqAM8BD+QejbvPkIDdnR9Ft5/A7bffx9KlG5k587WM56O6basrRASw7wfb1c9ep425JGNfsqrAW4BtxJB5AwbrjMVjnbNnxb7bEUkCjvl6VQUJImC7J511VlvOO+9sXnnFblipiIAIhCTlcQGni9a2Y7H3fk3ERwO2XeGnwKPh/YLdku0ygE0gqwaUCAeV3WcR45WAIyLSAS4isG3bdho0uIB27e6mZ8/rXRSZQhGB5BHwg4CTRy/3LRcKz8geA1wYPj27zyLWLAFHRKQDXEbg11//4IwzGjFkyKvUrduIFSugcmUoX95lgSocEUgQAQk4QaCzNGPPrtMFnP6l7D47ZHQScHISp1ZjIzBv3g9cfPEINm9+iiJF8pGSAoMHQyf/r9kRGzid7UsCEnDktA4EWgEnhZ8nf5vplCrAcKAssCm8QcSSyFWGJo9JwFGA0iH+IrB2LVSokMrevfsXr7C9G1at0kjYX5lWb6IhIAFHptQwvO60rap1OZBZwJ8DI4A3gLZAN6BOeBb24CxVT8r0upQEHJm7jvAhgWnToGVL2L59f+cOOww++QQaN/Zhh9UlEciBgAQc/eWxJouAbbKXvWdsk6hs1yWbdb0OOC/8eU4150XA9wH2J1RKly59/MaNNtdMRQS8Q8BGwBUrwp5M+5QVLLiPVavy6Vmwd9KoSB0iIAFHDzKrgM8OvwJlM5rTy9dAd8BGxocqJl97NWohcC+wOPw+c9bPcoxMz4CjT5yOdBeBIUOgSxcoVGgfO3fupWHD95k+3V5kCGbxy7rYwcxebL2WgKPn55SAo28xhyMlYEcwqpIkEbCR8MqV9tw3hSuuuJQLLmjN88/bzqHBKrYudsWBFUnZm6LtDIOV+lBvJeDok+7kLejoWz3EkRJwzAhVgUsI/PPPRqpVa8TgwU9w1VU23zE4RetiByfX2fVUAo4+/1kFbGdOD0/ASp+EZbefa0dfZd6PlIDzzk5nuo/Azz//TN26TRkx4m2aNUvIt5ArIGhdbFekIWlBeF3AJwI2d9JmIservAq0CO8r/Hd4ren0/YXt+a/J19aj3gzcDCyKVyCZ65WAE0FZbSSSwGefzadNm6tDGzdUqmRv/QWjaF3sYOTZjyPgKwFbVapA0FIoAQct48Ho75gxE3j44YeZO3cmZcocHoxOq5eBJeD1EbAEHNhLVx33K4H77hvEuHFj+fHHiRlbGPq1r+pXsAm4VcCrokxLMeAojYCjpKXDRMAjBK67risFC27ljTf+L2MLQ4+ErjBFIGoCbhVwSvhZqr1Xm1OpGF7SUbego065DhQB9xOwLQzPOKM1F1zQhMGDM9afcX/gilAEckHArQKeD6wOL++YU3d0CzoXydahIuAlAps2baZ69YYMGvQEl112mZdCd32s27dv56mnnuKjjz6iePHi2C88Rx55JC1atGDZsmXcdtttnHHGGa7px5tvvsnTTz/N0qVLM2I67rjjaNOmDV27dmX8+PEMHjyYFbbFFnDYYYdRvnz50COMTZs2UaZMGWrUqEGrVq1C5+TLd+jt63/99VdefPFFPv30UxYutPWS4lfcKuBXgEsBm+UcScDvAvtXdo8fK1fVrElYrkqHgokTgTVr1lCzZlMmTBjHeeedGadWglXtzp07adiwISVKlODdd98NidfKV199xc033xwS8KJFiw4QsEnuueeeSzqoe+65JyTH008/nQULFhwwR2D06NFce+21nHPOOXz55Zfkz5+mhZSUlJBMH330Ub799lvq1avHiBEjqGx7YWYp//d//xeS+ccff8xJJ52EXX/xLG4VcCXgdGB8hM4fBtiazD/HE5Ib65aA3ZgVxRQPAhMmzKFjx5tZvHgmxxxzTDyaCFSdL730EnfffTdff/01derY3jH7y2+//RYS76xZszIEPHXq1NCIeJVtWZXkMmzYMG655Rauvvpq3n777QOi+fzzz2natCkNGjRg9mzbO+fAsmvXLtq3b8+4ceOoWLEic+fOpWxZ28juwJKamsoRRxwR+hNUASc5ze5vXgJ2f44UoXMERo4cxWOPvcTixZ9TrJj93q2SVwJ33HEHr776Km+99VZISFnLf/7zHzp16hQS8Pz580NSO/zww+Muo2j688Ybb4RG6ddddx0jR4484JTp06dzwQUXHFLAdrDderd+rV69mquuuooxY+wt1oPLiSeeGBpBS8DRZCWAx0jAAUx6wLt89dWPsnr1SubOHaWZ0TFcC88//3zouak9+7XnpjfeeOMBtY0dO5aqVauGZGW3bWfMmBH6et26dUN/m7yrVasWulVtX//zzz9Dz1m3bNnC9ddfH/qsQIECbN26NXQr1+r75JNP2LBhA126dAnd9m7evDnvvfdernsRq4CtQbuFbbeyTbC//PILxx9//EFxnHzyyaHPJOAD0VwBTAa25TpzPjtBAvZZQtWdiATs1uAVV1zNsceeycsvPxrxeB2QPYFt27ZRv359vv/++9AB559/Pn369AmNHLMr2cnIZGojSXuWnC7SBx98kP79+zN06NDQbeJ169aFRpotW7bk33//pVu3bqHR9AsvvBAS9IcffpjrFDkhYPvF4dRTTw21baNoG01nLRJw9qnZCxi55bnOnM9OkIB9llB1JyoC27Ztp0KF8+nd+3/cdlvLqM7RQQcTsL3EbTRqE5fSi41KTcQ2WzhzyU5GNqq98sor6d27d2jlMit2u7p27dp07tw5NLJOLyb2OXPmhCZN1axpu67mvTghYPsFxCagWbH+du9uS/gfWCTg7HOUCpwiAWO3TfbZdHkVEQgagR9+WM1FF13IxIkTOeOMqkHrvqP9/eKLL0ICSp+0ZLdl7fZ03759Q6NUK9nJyG43263cDh06hG7h2qYSNrnJpHzTTTfx+uuvZ8TZuHHj0G1sOybW4oSAd+zYQbFitoYToRH7/fffLwFHmRgJOAxKAo7yitFhviTwzjtT6NLlflav/oJSpUr6so+J7NTkyZPp2bNnaJRq5aGHHgpJ+FACTo/NbjO/8sorodeWzjzzTP73v/+FpGyiTC9uE/BPP/0UesZtxW6f2y8NWYtGwNlffRKwBJzIn0tqy8UEOnXqx4oVs/n00/cpVKiQiyN1V2i2oMWdd95JyZIH/uJiz9jtOa2NCm0hC3tuW6RIkWxHwNajAQMGMHDgQF5++WUuvfTS0ISlChUquF7A9guD3Sa3vtkvEPa6kQQc3TUqAUvA0V0pOsr3BOyW5g033MOPP/7DV1+NyLhl6vuOx9hBe25bvXp1WrdufVBNxtTekTWZrl+/PvTedXajQXvGaxKfMGFCaJKVFS8I2BblsL4vX7489C60/QKRXdEIWCPgHL/NdAs6xp9COt0XBGzU1rDhrRx9dAHef18bN0STVBOwvR5ki20ULFjwoFNshSh7FmySslKpUiV2794demUnvdiI157B28SrWrVqhT62Gc8mb3utafjw4RnHpt+C3rNnT8y/JNkMa3tH2d5ftveYM5dp06bRpEmTQ74HbMttduzYMbQKlt0ut4lh9ipWdsXeA7blKn/+Ob5rPLl1JaxDXUcaAWsEHM3PGB0TIAL2g7V16+tITS3Hxx8/p3eEI+TeBGzv6l500UUMGjQoY0lGGx3+97//Db0mNGnSpNArRlbs73nz5vHXX3+F/titaVub+dlnnw0tfGEzie3VHluZys6z1bXslSSTsYnOXvlZuXIlM2fODL3yFEtJX4rSpG8rWWX+BeKdd97hmmuuyVbAtsymTTazyWCXX3556Bl16dKlsw3F3n8++uijQ7+E2C3qQ0k6ln6knysBO0ExCXVoBJwE6GrStQRshHbmmVdwwQV1GTw47bUYlewJpAs4/au2qIY9B7XFNM4666yQnE877bSMk21kecMNN4RmO9vo1l5fsteYbEUqW6bSNj2w937tlrQt1mELc5iUTcC2MEf6yLlo0aKh820hj9wWE36/fv1YvHhxxqm2jKQtSWmztu2dYlti00bhVkygtpazPee1WO1vex3K4sm6/GbmWKxvVtfvv/8e+tj6bCtmxWsdbAk4t1eCS46XgF2SCIXhGgJbt26jXr3G3HJLd+699+CZra4JVIGIQJiA1wTcARgLbAl6BiXgoF8B6n92BFau/JWaNRszbtw7NG16tiCJgKsJeE3AroaZyOAk4ETSVlteIjBlyjw6dbqW6dNncPLJx3kpdMUaMAISsEcTLgF7NHEKOyEEBgwYwzPP9GflyukUL5626pGKCLiNgFcFXNj2WXYbzETGIwEnkrba8iKBm2/+Hxs3LuH999/O2Jzdi/3wY8y29OUjjzwSddfSd2CK+gSPHOhVAfcGugD2JvmsLKzLAbaI6W8eyUGewpSA84RNJwWIgC0qUb9+eypUOIW33no8QD1XV71CwKsC/gI4HTgG2JUJdn+ga/j/84AbgJ+8kozcxCkB54aWjg0qge3bd1CtWmP693+Qq69uG1QM6rdLCXhVwPaS1gp7RzwT1zMA2+DSttz4BzgSsGVMqgNbXco/z2FJwHlGpxMDRmDt2rWcfXZjJk36hFq1qgWs9+qumwl4VcA26rXXkdpngvsUYBs73gvYAp/2Nv4TQE/gaTcnIS+xScB5oaZzgkpgxIhPeeCBh1izZk7GVnRBZaF+u4eAVwW8CZgGXJ4J5XfACeHb0nuA/OERsD0Lruse5M5EIgE7w1G1BIfAI4/8jwULVvDxxyO0XGVw0u7qnnpVwPOBMkCFMN3yYdmOAa7JRPxDoBFwuKuzkIfgJOA8QNMpgSZgGzfUrt2eMmXKMmXKS5JwoK8Gd3TeqwJOv7V8BzAEeAXoBNwKvJ4J7SjAZl4UcQdu56KQgJ1jqZqCQ8DWjG7Z8trQNnvDh78oCQcn9a7sqVcFbKPfpcBRwGagFLANOAn4NxPpGUAVwHfL4UjArvx+UlAeIGASrlPHJHwskya9IAl7IGd+DdGrArZ82GtINvq157vbgdsBG/GmF1us4++wqM/1WwIlYL9lVP1JJAGTcKNGl9OiRWMefvjBRDattkQgg4CXBZzeCVtnzmZF782S1zbA+8ALmd4N9k3qJWDfpFIdSRKBLVu2ULlyIx5++BHuueeKJEWhZoNMwA8CPlT+mgJXAW8D0/2WZAnYbxlVf5JBYMWKtTRocAEffTQ6x31ikxFbvNocPRr69k2rfc0aKFECypZN+3+3btA+88udDgTxySdgq06mpsKePfDgg9DB9rUDJk1K+1pKChQrBrZV8Flnwc6dcM01sHQpHHaY7e8LL78MlSunnXfRRbB+PeTPDyVLwsCBULPmwcH+/Tc0NROEy/btsGoVbNgAZcpEX48DGLKtws8CjhczV9QrAbsiDQrCBwS++WY+zZpdzdy5n1O16ok+6FH0XTjnHOjRA9rY/cI4lH374MgjYfp0OPPMNOGfcgr8+WeajE2oM2fC6afDrFnQuTMsXpwm4M8/h0sugXz54KWX4L330uqxsnEjHB5+t2XcOHj8cfjOXkSNUPr3hxkzYMKE6Ov59lto3RquvTaNlcm7SRPo0wfatoUiMUzxlYAjZcylX5eAXZoYheVJAoMGfUCfPv9j6dIZlCplczr9X0yANnr84Qc4+eT9/TW5mXTsM5NcuujyQsQEbKNrk2TDhvD992lSXb067d8mteXL99ds6E2ytWod2No336TJzgSetbzxBjz/fFrMkcqpp6aJ8/LMK0iET8qpnsaN00bqF1wATzwBN910ILNI7R7q6xJwXskl+TwJOMkJUPO+I/DMMwN4883PWLBgPAULFvRd/7J2aNEiaNQI/rGFezOVTz+Fd95Jk7Pdjj43mymsV18NP/6YPSIbXZa3lRnCZerUtNvJxYvDv//C2LFw4YWwaRNUrJg2Gq1fH8aPTxtpvv8+XJHlkfwNN6TdMn7BZvSEy403wjRbjgmw29zVbdHhHMqcOWn1/vorZE5vNPWYgLt2TRul29/HH+/M5SEBO8Mx4bVIwAlHrgZ9TsB2T2rSpDNFihRi0qQXfd5bGD4cRoyAzz7LvYCjhWOjbJNtr15pI+B586BVKzD528jYBGq3j7duhXr10v5vI1Q7Jr089VSapC1Oe06ctVg/7BcGk3BO5ZZb0m6H9+uX/VE51WMC3rEDjjsOxoyBQoWiJZDzcRKwMxwTXosEnHDkajAABOz1pIsuakGdOi3p1+8eX/f43nvTRPLMMwd2M5pb0NGOgO3WcdbbzHXqgEm1WbMD2921C8qVS5N0+mQre2b79ttgo+icboXbRC0b2Zpgsysm+GOPTavbnkEfqhyqHhNw9+5pv7AYHxul27GxFgk4VoJJOl8CThJ4Net7Av/+u4nKlc+nT5+nuO22lr7tr41I77gjTZDxKn/8kSbTr78Ge/66YgXYxC97XnviibBuXZoYrdgzVnsebXKz8uyzMGpUmnyPOGJ/hCZAm81so1ErH3wAd90Fa9emTdjKrrz2Grz+Osyenbd60p8B2+SrLl3SZmfbqLx06djIScCx8Uva2RJw0tCr4QAQWLZsDc2aNeXtt9+nQYMavuuxTY6yEeXcuTmPCJ3ouL32ZCNee2XIXkWymcTp0u/UKe25qt2qtlvQL76YFpeNZu05sj0jtmfRVmy2scX7889w1VVpt4StzqOOAhsp1win6dZb025hZ76Nbc+Yra2bb97fo0j1pB85f37apC17Ht6zZ5r8bUa33UK3XxratYPCtuxTHooEnAdobjhFAnZDFhSDnwlMnvwVV111A99/P4OTT/bdarZ+Tp1n+iYBeyZVBwYqAXs0cQrbUwT69x/Da6/145tvZlDcpvGqiICDBCRgB2EmsioJOJG01VaQCdx//1NMmjSX778fS4ECBYKMQn13mIAE7DDQRFUnASeKtNoJOgF7Pem88zpSrdqRDBvWP+g41H8HCXhVwNWAQ7wG7iAdF1clAbs4OQrNdwRSUlKoX/9imja9hr59beM1FRGInYBXBZwa3gd4AfBN+M88YHXsSLxRgwTsjTwpSv8Q+PPPf6hW7Xxefvk5rr76Iv90TD1JGgGvCngJYPtiZF0v7t9MQjYxzwE2JI1uHBuWgOMIV1WLwCEILF26gosvvpgxYz6kXr0zxEkEYiLgVQFbp58EugN/AasAW5ekKlAU2JeJii3R/QrwGmAjZ18UCdgXaVQnPEhg3LhZ3HTTrSxZMpMTTjjGgz1QyG4h4FUBdwZeAnoCtpBaulhtRNwa6A2kPye2kXJ+wG5X27I2fyQRfhlgSji2EuE4zgGeA2wNl89tQZho4pOAo6GkY0QgPgQGDBjJqFEvMXv2NIoVc2BNwviEqVpdTsCrAl4M7ARqH4KvjYJtx0fbV+zq8Ej5NmAhYHt77ElSXmwJb1vXZQxwYTgG+2x3+N+2LLrtzLk5UnwScCRC+roIxJfAtdf+l6VLf+DDD99m1ar8oSUXM+8CFN/WVbsfCHhVwDuAcUBOq5geFZ6U9RAwGHgeuBvoBAxLcvKmZhJweij2guEbgC2WFvEXBAk4yRlU84EnYK8nVakyhNWrb6Fo0QKkpMDgwWlLHqqIQDQEvCrgX4H1OYyA0/v+MWC3fesBhwO/AF8BuZnCOBCwzbFOAmoCmbd9rgIMB8oCm4CbAJsgFqlkFbD9IvE4MDn8S0Kk85GAIyLSASIQVwK2+H/FivvYs2f/DgC2z+yqVRoJxxW8jyr3qoBtpHhDWKRZdrM8IDu2r4bd6k3fs2Jm+PlrbmZONAxP8rJ9NC7PImB7ZjsiPHJtC3QD6gCnhUfdmYOZBDwd/iC7EbA9p343LOJFka4xCTgSIX1dBOJLwPaubdkybXH+9GJb1Nm+tLZ7jooIRCLgVQHbjo42ErVL/0bgo2w6agu3Lg8/Bw7vp8FbwJW2sUYkMNl8fU0WAR8NrAiPsO2Wsf0avA44L/x5Tk1kFrDFsit8sP1iYbO7f8rm5PsA+xMqpUuXPn6j7culIgIikBQCaSPgtJ180kuBAvtYvTqfngUnJSPea9SrAjbS7cMjT5v5bCPRkYAtxrE1/I7wo4CNXu1r6Vs/TwTqApl2l4w6aVkFfDaEhG6zrdPL1+EJX9bmoYrJ125l24Swe4FTgTvDM7VthK5Z0FGnRAeKQHIJDBmStj+sbWyfkpJK4cKP8csv91LW9qpTEYEIBLwsYOtag/D7vfb+b+Z3f+1rNiK1mdKNwmK2z0yi9pmNoHNbnBJwbtvN9njdgnYEoyoRgZgJ2Eh45UqoVAkmTnyNgQPfYP78qRSxDWxVRCAHAl4VsI0g7Tmp3fyxZ6f2qtEV4UlZdmvYXuOZBTwRPs4Q2IjVRsj/B9yRh6vCyVvQeWj+wFMk4JgRqgIRiAuBNm2689tvvzJ37pvky7d/glZcGlOlnibgVQHbwhs2m7l+Lunb5CsTtj2rzW3JKmA7f3r4Nrg9u7VJWLYy16HeTc5tezkeLwE7ilOViYBjBFJTU7n00nacfnp1Bgz4r2P1qiL/EfCqgG32kU28uj4BKXkVaAGUA/4GtoSfMVvT9vzX5HtkeNRt7/BGnMHsRMwSsBMUVYcIxIfA9u3bqVTpAu6++x569rwuPo2oVs8T8KqAZwB7gSaez0AeOyAB5xGcThOBBBFYu3Y99eo1YvToYZx/vk1XURGBAwl4VcA28rXVrGw7EnvVKHBFAg5cytVhDxL4/vvFnHdea7744lOqV6/kwR4o5HgS8KqAbWEN293IXuG5LLxQRjw5ua5uCdh1KVFAIpAtgTfemMwjj9zPokWzOOKIvLwBKbB+JeBVAdvt59+B48PPZAeENzhY5tdEZe2XBByUTKuffiAwcODLDBz4LkuXTqJw4cJ+6JL64AABrwrYJmHZTkfpJf0dYJvdPD/8utE3gP2x/YJ9VyRg36VUHfI5gRYt7mPXrk1MmTJUryf5PNfRds+rArb+2QMVe7e3Vqa/0+/vZF6UwzZgqBAtEK8cJwF7JVOKUwTSCOzdu5cWLdpw0kn1efVVe2NRJegEvCzg7HJ3cljG6WK2d3JNyrbVn6+KBOyrdKozASGwZctWKlVqRI8ePeja1ZYOUAkyAb8JOLtcnhjehtBXeZaAfZVOdSZABFat+o0mTRrzyisjad783AD1XF3NSiAIAvZl1iVgX6ZVnQoIgWnTFtKqVVvmzfuMU06xG3cqQSTgdQHb/ryXALb1yL/AB4fYmtB3uZWAfZdSdShgBF59dQLPPfcI7747i7/+KkXlymgbw4BdA14VsD3TfT/8DnD6auc28Wo40DGcQ5PyP4CtG+27IgH7LqXqUAAJXHzxFD799AKKFStASko+Bg+GTp0CCCKgXfaqgB8CngZWhzewt40SbJ9dW5c5XcD29wtA6/CewL5KsQTsq3SqMwEkYNsYVqy4jz179u+YVLAgrFqlkXBQLgevCngxYJOrbF9fW5DDio10MwvYVstaD7wCdPVbQiVgv2VU/QkagWnToGVL2L59f88POww++QQaNw4ajWD216sC3hHeCtCe/6aXrAK2z23LwoKJ2iIwkZeQBJxI2mpLBJwnkDYChj22q3m4FCiQyurV+fUs2HncrqzRqwK2lbA+B66IIOAxgG1DYktW+qpIwL5KpzoTUAJDhkCXLlCoEOzevY/ixXsyduxVNGli6wup+J2AVwX8JXAMUDGCgEcBVwJF/ZZICdhvGVV/gkrARsIrV0KlSrBkyde0a3cd3347jYoVTwgqksD026sC7gH0Bmy+oG1LaCW7W9DTw0tVZl432hfJlYB9kUZ1QgQOIjBs2Pv07/8kc+fOpGTJEiLkYwJeFbAtL2k7H9lEqw7AO9kI2H59XAEsCS9P6as0SsC+Sqc6IwIHELj11r588cVsFi/+gAIFfLeSrrIdJuBVAVv4DcOLbhQH5oSf9doz31uBM4BBQA3gCeBxv2VcAvZbRtUfEdhP7WSaYwAAIABJREFUYN++fbRs2YmKFUvy4ovPCY1PCXhZwJYSE+yIsHDt/5l3QbKX61YBNcN7BvsqhRKwr9KpzojAQQR2797Nqac25/LLr6R//y4i5EMCXhewpSQ/cE14spXJ9mhgE/AJ8Gj4XWDfpU4C9l1K1SEROIjA33//S/XqDRk27BmaN28uQj4j4AcB+ywl0XVHAo6Ok44SAa8TWLlyJXXqXMSECR/QoEF1r3dH8WciIAF79HKQgD2aOIUtAnkg8OGHX9Cx480sWTKTcuXK5aEGneJGAhKwG7MSRUwScBSQdIgI+IjAW2+N5tFHn+f776dRvHgxH/UsuF3xioB7At8BH8eQqhbAWcBTMdThmlMlYNekQoGIQMIIXH99LxYv/p4FC8aQP79Nf1HxMgGvCDi7RTZyy/114EbAFy/VScC5Tb+OFwHvE7DXk6655kaKFTue11+3DeFUvExAAvZo9iRgjyZOYYtAjAR27dpF5crN6NChA7173xJjbTo9mQS8JOCtwF8xwCoL2KIdGgHHAFGnioAIJJ/A+vV/Ubfu+QwaNIgWLZokPyBFkCcCXhJwnjqY5SRbqEMCdoKk6hABEUgqgUWLfqR+/UuZNu1jate2rdFVvEbAKwI+yUGwPztYV9Kq0i3opKFXwyLgGgLvvDONBx7owoIFMznqqKNcE5cCiY6AVwQcXW8CdJQEHKBkq6sikAOBQYNep2/f11i+fCpFi/pu51Vf514C9mh6JWCPJk5hi0AcCLRt25N///2ZqVNHki+fLYOv4gUCErAXspRNjBKwRxOnsEUgDgRSU1Np1epqDj/8dEaO9N3mb3Eg5o4qJWB35CHXUUjAuUamE0TA1wS2b99BpUoX0LXr3Tz00HW+7qtfOicBezSTErBHE6ewRSCOBH799Q8aNmzIs8++xuWXnxfHllS1EwQkYCcoJqEOCTgJ0NWkCHiAwJw5i7n44tbMmfMp1atX8kDEwQ1RAvZo7iVgjyZOYYtAAgiMGDGZJ5+8nzlzZnHkkUckoEU1kRcCEnBeqLngHAnYBUlQCCLgYgLdur3MmDHv8uOPkyhcuLCLIw1uaBKwR3MvAXs0cQpbBBJIoFWr+yhVahNvvjlUryclkHu0TUnA0ZJy2XESsMsSonBEwIUE9u7dy7nntqFGjQYMHdrNhREGO6QgCbg8UD38p6/X0y4Bez2Dil8EEkNgy5atVK3amGef7U779m0T06haiYqAHwVcOpNo04V7BlAKsCVidgGHRUXHxQdJwC5OjkITAZcR+O2336hduzGvvz6K5s3PcVl0wQ3HLwJuCDwInAUcnymdfwDfA9+F/7Z//wDs8XrKJWCvZ1Dxi0BiCUybtpDWra/iu+8+o0IFJ/e3SWw//NSaXwS8BigBjAV+Am4HbO793cBbfkpYel8kYD9mVX0SgfgSeO+98XTv/ijz5s3kiCPsZqFKMgn4RcCpQAtgYhimbQnSC7gXmATcAfyeTNBOty0BO01U9YlAMAjcdddzTJ48mR9++IiCBQsGo9Mu7aVfBFwX+BH4NwvnmsBQwJaDsSmAr7o0D7kOSwLONTKdIAIiAOzbt4/27e+kWDF47bVBej0piVeFXwScE8L8wH2AbREyD7gVWJlE5o40LQE7glGViEAgCezZs4dTT23BJZdcysCB/wkkAzd02isCtpGtvTr0dAzQKgKvAA2A4jHU44pTJWBXpEFBiIBnCfz77yZq1GjIiy/2plWryzzbDy8H7hUB2zPeN4CO2cC2V4v25SIJNwBv5uJ4Vx4qAbsyLQpKBDxFYM2aNdSocSHjxr3LBRfYEzuVRBLwg4A3h0fGTyUSXB7bKgNMAaqFZ21bNScDX4afYf8C3BhN3RJwNJR0jAiIQCQCEyfO5cYbr+fbb6dz/PGZ3+KMdKa+HisBPwg4p9FxrHycPr8QUBIYA1wYrtwE/Ej42XTU7UnAUaPSgSIgAhEI2KYN99/fl6VLp1OypL3RqZIIAhJwIigf3MbULAKeBdi7zPaMelQ0IUnA0VDSMSIgAtES6NixD/Pnf8WCBWMpUKBAtKfpuBgISMCR4Q0EWgG2dIw9JPk20ylVgOFAWWATcBOwJHKVZBZwEcCudnuOPRloA/wdqQ4JOBIhfV0ERCA3BOz1pBtv7MiuXWUYM2ZAbk7VsXkkIAFHBmfLXK4CZgOXZxHw58CI8AQxW+Xc3jWuA5wGDM5StS0Ikj6LO7OAMx/2JPAh8HWksCTgSIT0dREQgdwSSElJoWrVi2nX7mr69bP1i1TiScBLAl4cXkjD/rY/6aPERD0DtlvEmQV8NLACsIlVtra0zcZeB5wX/jynvGUWsD1w2Ro+3yRto2irJ2uxd5ntT6iULl36+I0bN8bz2lDdIiACASTw55//UL/++Tz55HO0a3dRAAkkrsteErBRyfy60YawiJuGl5u0BTbiudxkVgGfHV5n2mY0pxcbuXYHbGR8qGLytVvZC8NLZR4H2Azu3cB7QFT3fjQCTtw3iVoSgaARWLRoBfXrX8ynn35IvXq2mZxKPAh4RcB2e9eEZ39qhUedmXmki9l2P/omvOJV+t9/OQTOKQE7Eo4E7AhGVSICInAIAuPGzeLee2/lq69mcuyxx4hTHAh4RcBZu24TotKFnC5lmwiVXjKPlO3d2goOsHPyFnTM4UjAMSNUBSIgAhEIvPjiSPr0GcSKFZ9TrJjnt1F3Xb69KuDsQJbPMko2MdtzWpOxE3PqswrYYpgenoBlq3TZKN1uP9dORJYl4ERQVhsiIALt2/+Xv/76gcmT3yZ/fltaX8UpAn4ScHZMbFkXE/H4GIDZDkq21WG58MSvLUDlcH32/NfkeyRgK3LdDCyKoa2oT5WAo0alA0VABGIgYK8ntWx5HcWKVeDdd+1FDRWnCPhdwE5xcl09ErDrUqKARMC3BHbs2EmVKk25995OPPCAvaih4gQBCdgJikmoQwJOAnQ1KQIBJrB+/Z/Uq9eQfv1e5qqrGgeYhHNdl4CdY5nQmiTghOJWYyIgAvZ6ybwfaNKkJTNnTqRmzapiEiMBPwrYNjw4PfyubQ3gLMB3v65JwDFe+TpdBEQgTwTeffczevS4my++mMkxx2R++SRP1QX6JK8LuDRgkrU/triF/X0KYBK2kr5XsBOzoF11oUjArkqHghGBQBF4+OGhvPnmCH76aQpFithy9ip5IeAlAdtM4waZhGuytfeB04vJNr3YGo220tQCYD7wdl7guPkcCdjN2VFsIuB/Atdc0438+dcxatRw8uXL/OPX/313qodeEXB74PUsI9t0BrYL0dywbE249sc2T/B1kYB9nV51TgRcTyA1NZU6da7ijDNqMny4bWmuklsCXhHwD4C9c2ubFMwMj26bh5/tmoDvBkbmtvNePl4C9nL2FLsI+IPAtm3bOfXUxvTp05XrrrNxkkpuCHhFwCnhjRfqAvbv9NIVsDfD7SHEu4DtnxWILYIk4Nxc5jpWBEQgXgTWrVtHjRqNGDr0DS67rH68mvFlvV4R8Pbwc9yO2WThVODN8CYNv4VXo7Idh3xdJGBfp1edEwFPEZg9+3tatGjDwoVTqFixoqdiT2awXhGwLUBqk7D+PASsgsBj4bWY7dhBwEPArmTCjWfbEnA86apuERCB3BL4+ONP6NKlGwsXzqJMmcNze3ogj/eKgKNNzjnACMDeEF8K3BB+Xhzt+Z45TgL2TKoUqAgEhkDXri8xfvyHLFv2CYUKpb8NGpju57qjfhOwASgKPAN0AfYA/wtveJ9rOG4+QQJ2c3YUmwgEl8Ctt/6HzZu38847/6fXkyJcBn4UcHqXLwy/unScQ9sRuuo7SgJ2VToUjAiIQJjA3r17Of301jRu3JhXXnlAXHIg4GcBW7ftQcSL4VvRvroQJGBfpVOdEQFfEdi8eQu1ajWkT5/HuOqqNr7qm5Od8buAnWTlqrokYFelQ8GIgAhkIfDLL2upXr0J7703mmbNaotPNgQkYI9eFhKwRxOnsEUgQASmTp3PddddzTffTKN8+fIB6nl0XZWAo+PkuqMkYNelRAGJgAhkQ+D99z/gzjv/x/LlMylVqqQYZSIgAXv0cpCAPZo4hS0CASRwxx39+eqraXzzzYcULGjLNqgYAQnYo9eBBOzRxClsEQgggX379tGx4+38+WdRPvpoYAAJZN9lCdijl4IE7NHEKWwRCCiB3bt3c8opl9KqVSuee872z1GRgD16DUjAHk2cwhaBABP455+N1K9/Pt279+Wmmy4NMIm0rkvAHr0EJGCPJk5hi0DACSxduppzz72QqVMncO65pwWahgTs0fRLwB5NnMIWARHgpZfeZ9SoF5kx41MKFy4cWCISsEdTLwF7NHEKWwREIETgiiu68scfG5k9e1hg14yWgD36zSABezRxClsERCBEwNaMrlnzSi6/vA69ej0cSCoSsEfTLgF7NHEKWwREIIPAtm3bqFGjCbfeei/durUPHBkJ2KMpl4A9mjiFLQIicACBtWvXU716I8aOHUqTJucHio4E7NF0S8AeTZzCFgEROIjA4sVLaNasFRMnTqJGjSqBISQBezTVErBHE6ewRUAEsiUwatRU7rrrHn78cSZHH102EJQkYI+mWQL2aOIUtgiIwCEJPP/8a4wY8Tpz5kylaNGiviclAXs0xRKwRxOnsEVABHIk0LZtT1avXs28eaPInz+/r2lJwB5NrwTs0cQpbBEQgRwJpKam0qzZtdSuXZm+fXv7mpYE7NH0SsAeTZzCFgERiEhg586dVK16ITfd1JFevTpGPN6rB0jAHs2cBOzRxClsERCBqAisX/8XtWo15M03X6Rp06ZRneO1gyRgr2UsHK8E7NHEKWwREIGoCSxfvpxzz72ESZP8uXGDBBz1peCuAyVgd+VD0YiACMSHwNixs7j99ltZtGgG5cqVi08jSapVAk4S+FiblYBjJajzRUAEvEJg5MhRPProQBYvnkbx4sW8EnbEOCXgiIjceYAE7M68KCoREIH4ELjxxif47ruFLFjwLgUKFIhPIwmuVQJOMHCnmpOAnSKpekRABLxAYN++fVx//c2kph7J6NEDvBByxBgl4IiI3HmABOzOvCgqERCB+BFISUmhatXmXHllWwYM6BK/hhJUswScINBONyMBO01U9YmACHiBwN9//0vdug3p3bsvV199qRdCPmSMErBH0ycBezRxClsERCBmAsuWraZOnQuZMOF9GjeuEXN9yapAAk4W+RjblYBjBKjTRUAEPE1g4sS53Hrr9Xz11TTKlz/Bk32RgD2ZNpCAPZo4hS0CIuAYgREj3ueBB55kxYoZlCpV0rF6E1WRBJwo0g63IwE7DFTViYAIeJLAHXf0Z8GCz5kzZzwFCxb0VB8kYE+la3+wErBHE6ewRUAEHCVgryfdeGNn1q3Lz5Qpg8iXL5+j9cezMgk4nnTjWLcEHEe4qloERMBTBPbs2cNpp11Gy5bNePbZ+zwTuwTsmVQdGKgE7NHEKWwREIG4ENi0aTP16jWic+fHuPvuNnFpw+lKJWCniSaoPgk4QaDVjAiIgGcIrFz5KzVrXsCYMaNo3vwc18ctASc2RWWAKUA1oES46eZA9/C/qwO28eW3kcKSgCMR0tdFQASCSGD69IXceGNbPv30M0455WRXI5CAE5ueQoDNlR8DXJilaZs58DVgv7btixSWBByJkL4uAiIQVAJDhnxM9+49+OmnmZQpc7hrMUjAyUnN1GwEbOK9GegcTUgScDSUdIwIiEBQCTz44CC++mosn302kcKFC7sSgwQcOS0DgVbASUDNLLeHqwDDgbLAJuAmYEnkKslOwE8Cs4BJUZyvhTiigaRjREAEAk3giiu68scfm5g9+zVXvp4kAUe+PBsCq4DZwOVZBPw5MAJ4A2gLdAPqAKcBg7NUbWJ9OvxZdgKeAzQGUiKHpJWwomGkY0RABIJNYO/evdSqdSWtW9ehV6+HXQdDAo4+JWuyCPhoYAVgE6v2APYMdx1wXvjznGrOKuCKYTm3y+Eke7kt4wW30qVLH79x48boo9eRIiACIhBAAtu2beOssy6gU6eudOvW3lUEJODo05FVwGcDb4VnNKfXYpOobEazjYwPVUy+dit7IXAvsBjoCmwARkUbjp4BR0tKx4mACASdwNq166levRHjxr3GBRfYGMkdRQKOPg9OCTj6FnM4UgJ2BKMqEQERCAiBxYuX0KxZKyZOnESNGjZ9J/lFAo4+B07ego6+1UMcKQHHjFAViIAIBIzAqFFTueuue1i+fCZHHWVzZ5NbJODo+WcVsJ05PTwBK30Slt1+rh19lXk/UgLOOzudKQIiEFwCL7zwGsOHv8GcOVMoWrRoUkFIwJHxvwq0AMoBfwNbgMrh02xFK5PvkcDm8Hu8iyJXGfsREnDsDFWDCIhAMAm0bduTNWvW8PXXI8mfP3/SIEjASUMfW8MScGz8dLYIiEBwCaSmptKs2bXUrl2Zvn17Jw2EBJw09LE1LAHHxk9ni4AIBJvAzp07qVKlKTfffAu9enVMCgwJOCnYY29UAo6doWoQAREINoH16/+iVq2GvPnmizRtavvgJLZIwInl7VhrErBjKFWRCIhAgAksX76cc8+9hEmTJnDuubaIYeKKBJw41o62JAE7ilOViYAIBJjA++/P5I47OrFo0QzKlbP5tokpEnBiODveigTsOFJVKAIiEGACI0eO4tFHB7J48TSKFy+WEBIScEIwO9+IBOw8U9UoAiIQbAI33vgE3323kAUL3qVAgQJxhyEBxx1xfBqQgOPDVbWKgAgEl8C+ffu4/vqbSU09ktGjB8QdhAQcd8TxaUACjg9X1SoCIhBsAikpKVSt2py2bdvSv3+XuMKQgOOKN36VS8DxY6uaRUAEgk3g77//pW7dhjzxRF+uuebSuMGQgOOGNr4VS8Dx5avaRUAEgk1g2bLV1KlzIRMmvE/jxjXiAkMCjgvW+FcqAcefsVoQAREINoGJE+fSqdP1fPnlNMqXP8FxGBKw40gTU6EEnBjOakUERCDYBEaMeJ8HHniSFStmUKpUSUdhSMCO4kxcZRJw4lirJREQgWATuOOO/ixc+DlffDGeggULOgZDAnYMZWIrkoATy1utiYAIBJeAvZ50442dWbcuP1OmDCJfvnyOwJCAHcGY+Eok4MQzV4siIALBJbBnzx5OO+0yLrusGQMG3OcICAnYEYyJr0QCTjxztSgCIhBsAps2baZevUZ07vwYd9/dJmYYEnDMCJNTgQScHO5qVQREINgEVqxYS61aTRg/fhSNG58TEwwJOCZ8yTtZAk4ee7UsAiIQbAJz5y6kSJF81KgR2/vBErBHryMJ2KOJU9giIAIiECYgAXv0UpCAPZo4hS0CIiACErC3rwEJ2Nv5U/QiIAIioBGwR68BCdijiVPYIiACIqARsLevAQnY2/lT9CIgAiKgEbBHrwEJ2KOJU9giIAIioBGwt68BCdjb+VP0IiACIqARsEevAQnYo4lT2CIgAiKgEbC3rwEJ2Nv5U/QiIAIioBGwR68BCdijiVPYIiACIqARsLevAQnY2/lT9CIgAiKgEbBHrwEJ2KOJU9giIAIioBGwt68BCdjb+VP0IiACIqARsEevAQnYo4lT2CIgAiKgEbC3rwEJ2Nv5U/QiIAIioBGwR68BCdijiVPYIiACIqARsOevgV3AnzH24jBgRx7qiPa8nI471Nei/TzrcVn/XwLYmoe+RXNKtP3PWle054nbgeTiyc1ayq7+SNdX0K83ccv+usl85Ubzs+wooEg0P3QiHZMv0gH6uusIjAHa5SGqaM/L6bhDfS3az7Mel/X/vwIn5KFv0ZwSbf+z1hXteeJ2ILl4crOWsqs/0vUV9OtN3LK/bjJfudH+LIvmZ07EYyTgiIhcd8BVwLt5iCra83I67lBfi/bzrMdl/X88BRxt/7OijfY8cTuQXDy5WUvZ1R/p+gr69SZu2V83ma/caH+W5eFH8MGnSMCOYFQlDhGIp4AdCtGV1Yhb3tIibuKWNwIOnSUBOwRS1ThC4D7gWUdqClYl4pa3fIubuOWNgENnScAOgVQ1IiACIiACIpAbAhJwbmjpWBEQAREQARFwiIAE7BBIVSMCIiACIiACuSEgAeeGlo6NB4GBQCvgJKAm8G08GvFJnYdiVQUYDpQFNgE3AUt80ue8diMvrILI0WlOQWSY12sUCTjP6HSiQwQaAquA2cDlEnCOVA/F6nNgBPAG0BboBtRxKD9erSYvrILI0WlOQWSY5+8RCTjP6HSiwwTWSMBRE83M6mhgBVAG2AOhX6rXAeeFP4+6Up8eGC2rzQHn6ASnoDPM9beQBJxrZDohTgQk4OjBZmZ1NvAWUC3T6V8D3QEbjQS9RMvKbt0HmaMTnILOMNffaxJwrpHphDgRkICjBxvtD0sJGKJlFXR5OMEp6Ayj/w4OHykB5xqZTogTAQk4erDR3i60W9NBL9GyCvrtUyc4BZ1hrr/XJOBcI9MJcSIgAUcPNiur6eEJWOmTsOz2c+3oq/P1kblhFWSOTnEKMsNcfyNJwLlGphMcJvAq0AIoB/wNbAEqO9yGX6o7FCt7/mvyPRKwUcjNwCK/dDqP/cgLqyBydJpTEBnm8RJNmzGpIgIiIAIiIAIikGACEnCCgas5ERABERABETACErCuAxEQAREQARFIAgEJOAnQ1aQIiIAIiIAISMC6BkRABERABEQgCQQk4CRAV5MiIAIiIAIiIAHrGhABERABERCBJBCQgJMAXU2KgAiIgAiIgASsa0AE/EOgQHghk0LASqBXeIOB/2/vPEMlKaIwehBBxRxARTBiQjFiAhVzQASzGDCA/jAhBhQVRVFMmMDwR0RMYM4IZsyYBXPAgBEVcxZUPqiGoZnHm57dN9Xrnvqz87ar6tacqn3f1q1bt8f5hisCH7Uang6cM05nE2rzTuulFE8AW03ItmYk0JmAAtwZmQ0k0FsCeTXhhcB6wLrA38BywDdjjLgR4F+B20v7O4D7xuhrUk3OA5YtWdV2BBTgSZHXzlgEFOCxsNlIAr0mMF/ZvUaMdgfuHmO0jQB/AuTznFSy631cAZ6TpmzuHKsCPHfOu9/6/0/gZmBf4CzgzDG+rgI8BjSbSKALAQW4Cy3rSmDOIXAycD5wD7DbGMMeVYCXB7JL/hTI52ElL9lYorws4rtSIb97vgcWLc+2B44E1gf+KbvX48pOfh7g8PKSibWAX4AHgBNLH22b7oDHmHCbTJ6AAjx55lqUwEwTyL/rx0oA0rgu5FEFeM9yRnwXsMeQL7ZyCQhLQFc+NyVvvHof+Ly4i9NPzmx/A7YBFgPeAzYB7gXWKc/TfgdgfuBOIO0U4JleUfY/IwQU4BnBaqcSqEogO8dLBkaQ3Wd2m13KqAKcwKe8f/g04NwhBvYBbikivffA87jH4yZPeRbIsy/KzxHqN4AFiki/CRwC/FiebwQ8D/wLLAT83rLrDrjLTFu3GgEFuBp6DUtgRgisCbwC/FUEKi7e7CgTlNSljCrAjwDbAjsBDw4xcAFwEnBKcYk3VZq//6y4nb9ttX0BiNC+BmwG/NF6HrFOkNkqwIcKcJeptW5fCCjAfZkJxyGBWSeQ+7/PARsCBwMHFHft8cClHbsfRYDz+yNnunEXLwXkrLddHi3/AYjb+OGBh41wn9DarTdV3gbWKC7muJoHS+z+DCwILFI+Dz53B9xxsq1eh4ACXIe7ViUwEwSSeCPJMnJvNy7dxj18A3BQR4OjCPCq5Zz2Y2ClIf0PCvSSRaybahHuxcsVp5xTD5YIa9zNucecOu3d7+pAkm60z5WbPhTgjpNt9ToEFOA63LUqgdlNYGPgmZJ0Y+0idhHhW4HXSxBTF5ujCPB+JdNWEnTsNaTzJtCqLdAR67iNvypu5HbTzYGnyjnvpkP63R+4aci5sgLcZYatW52AAlx9ChyABGaZQIKVXgVWA3YeOIvN+egHZSe5MPBnB0ujCHACvRLwNVWKymOBy4C2QEesbwPuB3YdMqam3RXAMUOeXwTEdZ3gr5wlt4s74A4TbdV6BBTgeuy1LIHZReBy4GjgyvLnYL+Jfs4Zbc6FE5w1ahlFgJ8Etij3cS9udZxsXC8B2Y2fWtzhTZXGNZ4EIUkU0i5xmR9YIp+vG/I8AWUR2dwdzlmyAjzqrFqvVwQU4F5Nh4ORQGcC2wEPAe8CGwy5kpP7wFsDhwHXdOh9OgFOcoyfSiDUi0UQc4c3JdeeriqZuPJz7gfnnnBTEi2doKxdSkKN9rCaAKwk3Xir9bCdwGPY9Sp3wB0m2qr1CCjA9dhrWQKzSiA725zvLlOu6mTH2S6Nu3Yqd+5UY5hOgCOOuav7JRD3dgKlYj+fkzQj57sJoEp0dAT1+oFrSLlylKCspYGvWwNI+x9KQo5coUpWrMHSnCvnDDku9mFFAZ7VlWX7iRBQgCeCWSMSmBECN5arRmcAZ09hoQmUerq4i0cdyHQCnMQY15ZkGleXtzDF3Zzo5Zztxu18aDkfTprKI8pd5BWABGVNlbpyy5LxKkFY+dwuTQKPBJflswI86oxar3cEFODeTYkDkkAvCEwnwNlRH1WCoQazbvVh8O6A+zALjmFaAgrwtIisIIG5ksB07wNOKshcfcouNbvVPhTfB9yHWXAMIxNQgEdGZUUJzFUEGgEe/NLNdaNk3EomqnlLJqom+Ko2oCTnSJKOpuTlDtkNWyTQSwIKcC+nxUFJoNcEEm398pgJPnr9xRycBCZJQAGeJG1tSUACEpCABAoBBdilIAEJSEACEqhAQAGuAF2TEpCABCQgAQXYNSABCUhAAhKoQEABrgBdkxKQgAQkIAEF2DUgAQlIQAISqEBAAa4AXZMSkIAEJCABBdg1IAEJSEACEqhAQAGuAF2TEpCABCQgAQXYNSABCUhAAhKoQEABrgBdkxKQgAQkIAEF2DUgAQlIQAISqEBLUc1OAAAAzklEQVRAAa4AXZMSkIAEJCABBdg1IAEJSEACEqhAQAGuAF2TEpCABCQgAQXYNSABCUhAAhKoQEABrgBdkxKQgAQkIAEF2DUgAQlIQAISqEBAAa4AXZMSkIAEJCABBdg1IAEJSEACEqhAQAGuAF2TEpCABCQgAQXYNSABCUhAAhKoQEABrgBdkxKQgAQkIAEF2DUgAQlIQAISqEBAAa4AXZMSkIAEJCABBdg1IAEJSEACEqhAQAGuAF2TEpCABCQgAQXYNSABCUhAAhKoQOA/8m5MsjcCpb4AAAAASUVORK5CYII=\" width=\"599.9999910593034\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.ticker as mtick\n",
    "\n",
    "# SET UP AXIS\n",
    "ax = plt.figure(figsize=(6,6)).add_subplot(111)\n",
    "\n",
    "# PLOT PHOTOSPHERIC AND LITERATURE DATA\n",
    "tauA.plot_photosphere(ax, pointsize=4,lcolor='blue',pcolor='blue',lw=.5)\n",
    "tauA.plot_observedData(ax, pointsize=3,lw=.8)\n",
    "\n",
    "# ADDITIONAL PLOT STUFF\n",
    "plt.loglog()\n",
    "#ax.set_xlim(.2,100)\n",
    "#ax.set_ylim(1e-13,1e-6)\n",
    "ax.set_ylabel(r'$\\lambda F_{\\lambda}\\ [erg\\,\\ s^{-1}\\ cm^{-2}] $',fontsize=20, family='sans-serif')\n",
    "ax.set_xlabel(r'$\\lambda\\ [\\mu m]$', fontsize=20, family='sans-serif')\n",
    "plt.subplots_adjust(left=.19, right=.96, bottom=.15, top=.93)\n",
    "ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%d'))\n",
    "\n",
    "\n",
    "ax.annotate(r'{}'.format(tauA.starname), xy=(0.73, 0.53),xycoords='axes fraction',\n",
    "            fontsize=20, family='Times New Roman')\n",
    "\n",
    "ax.annotate(r'$T_* = {:.2f}$ K'.format(tauA.StarTemp), xy=(0.73, 0.47),\n",
    "            xycoords='axes fraction',fontsize=10, family='Times New Roman', color='blue')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 4. Calculate Excess Flux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HPACS100  excess flux  2.38746427586e-17 erg/s/cm2/ang\n",
      "HPACS160  excess flux  4.66817105066e-18 erg/s/cm2/ang\n",
      "W3  excess flux  2.14361397755e-16 erg/s/cm2/ang\n",
      "W4  excess flux  3.97925224579e-17 erg/s/cm2/ang\n"
     ]
    }
   ],
   "source": [
    "# There's a function that will calculate the excess flux at bands listed in mags4Dust\n",
    "tauA.calc_excessflux()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "/* Put everything inside the global mpl namespace */\n",
       "window.mpl = {};\n",
       "\n",
       "\n",
       "mpl.get_websocket_type = function() {\n",
       "    if (typeof(WebSocket) !== 'undefined') {\n",
       "        return WebSocket;\n",
       "    } else if (typeof(MozWebSocket) !== 'undefined') {\n",
       "        return MozWebSocket;\n",
       "    } else {\n",
       "        alert('Your browser does not have WebSocket support.' +\n",
       "              'Please try Chrome, Safari or Firefox ≥ 6. ' +\n",
       "              'Firefox 4 and 5 are also supported but you ' +\n",
       "              'have to enable WebSockets in about:config.');\n",
       "    };\n",
       "}\n",
       "\n",
       "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n",
       "    this.id = figure_id;\n",
       "\n",
       "    this.ws = websocket;\n",
       "\n",
       "    this.supports_binary = (this.ws.binaryType != undefined);\n",
       "\n",
       "    if (!this.supports_binary) {\n",
       "        var warnings = document.getElementById(\"mpl-warnings\");\n",
       "        if (warnings) {\n",
       "            warnings.style.display = 'block';\n",
       "            warnings.textContent = (\n",
       "                \"This browser does not support binary websocket messages. \" +\n",
       "                    \"Performance may be slow.\");\n",
       "        }\n",
       "    }\n",
       "\n",
       "    this.imageObj = new Image();\n",
       "\n",
       "    this.context = undefined;\n",
       "    this.message = undefined;\n",
       "    this.canvas = undefined;\n",
       "    this.rubberband_canvas = undefined;\n",
       "    this.rubberband_context = undefined;\n",
       "    this.format_dropdown = undefined;\n",
       "\n",
       "    this.image_mode = 'full';\n",
       "\n",
       "    this.root = $('<div/>');\n",
       "    this._root_extra_style(this.root)\n",
       "    this.root.attr('style', 'display: inline-block');\n",
       "\n",
       "    $(parent_element).append(this.root);\n",
       "\n",
       "    this._init_header(this);\n",
       "    this._init_canvas(this);\n",
       "    this._init_toolbar(this);\n",
       "\n",
       "    var fig = this;\n",
       "\n",
       "    this.waiting = false;\n",
       "\n",
       "    this.ws.onopen =  function () {\n",
       "            fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n",
       "            fig.send_message(\"send_image_mode\", {});\n",
       "            if (mpl.ratio != 1) {\n",
       "                fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n",
       "            }\n",
       "            fig.send_message(\"refresh\", {});\n",
       "        }\n",
       "\n",
       "    this.imageObj.onload = function() {\n",
       "            if (fig.image_mode == 'full') {\n",
       "                // Full images could contain transparency (where diff images\n",
       "                // almost always do), so we need to clear the canvas so that\n",
       "                // there is no ghosting.\n",
       "                fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n",
       "            }\n",
       "            fig.context.drawImage(fig.imageObj, 0, 0);\n",
       "        };\n",
       "\n",
       "    this.imageObj.onunload = function() {\n",
       "        fig.ws.close();\n",
       "    }\n",
       "\n",
       "    this.ws.onmessage = this._make_on_message_function(this);\n",
       "\n",
       "    this.ondownload = ondownload;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_header = function() {\n",
       "    var titlebar = $(\n",
       "        '<div class=\"ui-dialog-titlebar ui-widget-header ui-corner-all ' +\n",
       "        'ui-helper-clearfix\"/>');\n",
       "    var titletext = $(\n",
       "        '<div class=\"ui-dialog-title\" style=\"width: 100%; ' +\n",
       "        'text-align: center; padding: 3px;\"/>');\n",
       "    titlebar.append(titletext)\n",
       "    this.root.append(titlebar);\n",
       "    this.header = titletext[0];\n",
       "}\n",
       "\n",
       "\n",
       "\n",
       "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n",
       "\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n",
       "\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_canvas = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var canvas_div = $('<div/>');\n",
       "\n",
       "    canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n",
       "\n",
       "    function canvas_keyboard_event(event) {\n",
       "        return fig.key_event(event, event['data']);\n",
       "    }\n",
       "\n",
       "    canvas_div.keydown('key_press', canvas_keyboard_event);\n",
       "    canvas_div.keyup('key_release', canvas_keyboard_event);\n",
       "    this.canvas_div = canvas_div\n",
       "    this._canvas_extra_style(canvas_div)\n",
       "    this.root.append(canvas_div);\n",
       "\n",
       "    var canvas = $('<canvas/>');\n",
       "    canvas.addClass('mpl-canvas');\n",
       "    canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n",
       "\n",
       "    this.canvas = canvas[0];\n",
       "    this.context = canvas[0].getContext(\"2d\");\n",
       "\n",
       "    var backingStore = this.context.backingStorePixelRatio ||\n",
       "\tthis.context.webkitBackingStorePixelRatio ||\n",
       "\tthis.context.mozBackingStorePixelRatio ||\n",
       "\tthis.context.msBackingStorePixelRatio ||\n",
       "\tthis.context.oBackingStorePixelRatio ||\n",
       "\tthis.context.backingStorePixelRatio || 1;\n",
       "\n",
       "    mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n",
       "\n",
       "    var rubberband = $('<canvas/>');\n",
       "    rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n",
       "\n",
       "    var pass_mouse_events = true;\n",
       "\n",
       "    canvas_div.resizable({\n",
       "        start: function(event, ui) {\n",
       "            pass_mouse_events = false;\n",
       "        },\n",
       "        resize: function(event, ui) {\n",
       "            fig.request_resize(ui.size.width, ui.size.height);\n",
       "        },\n",
       "        stop: function(event, ui) {\n",
       "            pass_mouse_events = true;\n",
       "            fig.request_resize(ui.size.width, ui.size.height);\n",
       "        },\n",
       "    });\n",
       "\n",
       "    function mouse_event_fn(event) {\n",
       "        if (pass_mouse_events)\n",
       "            return fig.mouse_event(event, event['data']);\n",
       "    }\n",
       "\n",
       "    rubberband.mousedown('button_press', mouse_event_fn);\n",
       "    rubberband.mouseup('button_release', mouse_event_fn);\n",
       "    // Throttle sequential mouse events to 1 every 20ms.\n",
       "    rubberband.mousemove('motion_notify', mouse_event_fn);\n",
       "\n",
       "    rubberband.mouseenter('figure_enter', mouse_event_fn);\n",
       "    rubberband.mouseleave('figure_leave', mouse_event_fn);\n",
       "\n",
       "    canvas_div.on(\"wheel\", function (event) {\n",
       "        event = event.originalEvent;\n",
       "        event['data'] = 'scroll'\n",
       "        if (event.deltaY < 0) {\n",
       "            event.step = 1;\n",
       "        } else {\n",
       "            event.step = -1;\n",
       "        }\n",
       "        mouse_event_fn(event);\n",
       "    });\n",
       "\n",
       "    canvas_div.append(canvas);\n",
       "    canvas_div.append(rubberband);\n",
       "\n",
       "    this.rubberband = rubberband;\n",
       "    this.rubberband_canvas = rubberband[0];\n",
       "    this.rubberband_context = rubberband[0].getContext(\"2d\");\n",
       "    this.rubberband_context.strokeStyle = \"#000000\";\n",
       "\n",
       "    this._resize_canvas = function(width, height) {\n",
       "        // Keep the size of the canvas, canvas container, and rubber band\n",
       "        // canvas in synch.\n",
       "        canvas_div.css('width', width)\n",
       "        canvas_div.css('height', height)\n",
       "\n",
       "        canvas.attr('width', width * mpl.ratio);\n",
       "        canvas.attr('height', height * mpl.ratio);\n",
       "        canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n",
       "\n",
       "        rubberband.attr('width', width);\n",
       "        rubberband.attr('height', height);\n",
       "    }\n",
       "\n",
       "    // Set the figure to an initial 600x600px, this will subsequently be updated\n",
       "    // upon first draw.\n",
       "    this._resize_canvas(600, 600);\n",
       "\n",
       "    // Disable right mouse context menu.\n",
       "    $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n",
       "        return false;\n",
       "    });\n",
       "\n",
       "    function set_focus () {\n",
       "        canvas.focus();\n",
       "        canvas_div.focus();\n",
       "    }\n",
       "\n",
       "    window.setTimeout(set_focus, 100);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_toolbar = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var nav_element = $('<div/>')\n",
       "    nav_element.attr('style', 'width: 100%');\n",
       "    this.root.append(nav_element);\n",
       "\n",
       "    // Define a callback function for later on.\n",
       "    function toolbar_event(event) {\n",
       "        return fig.toolbar_button_onclick(event['data']);\n",
       "    }\n",
       "    function toolbar_mouse_event(event) {\n",
       "        return fig.toolbar_button_onmouseover(event['data']);\n",
       "    }\n",
       "\n",
       "    for(var toolbar_ind in mpl.toolbar_items) {\n",
       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
       "\n",
       "        if (!name) {\n",
       "            // put a spacer in here.\n",
       "            continue;\n",
       "        }\n",
       "        var button = $('<button/>');\n",
       "        button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' +\n",
       "                        'ui-button-icon-only');\n",
       "        button.attr('role', 'button');\n",
       "        button.attr('aria-disabled', 'false');\n",
       "        button.click(method_name, toolbar_event);\n",
       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
       "\n",
       "        var icon_img = $('<span/>');\n",
       "        icon_img.addClass('ui-button-icon-primary ui-icon');\n",
       "        icon_img.addClass(image);\n",
       "        icon_img.addClass('ui-corner-all');\n",
       "\n",
       "        var tooltip_span = $('<span/>');\n",
       "        tooltip_span.addClass('ui-button-text');\n",
       "        tooltip_span.html(tooltip);\n",
       "\n",
       "        button.append(icon_img);\n",
       "        button.append(tooltip_span);\n",
       "\n",
       "        nav_element.append(button);\n",
       "    }\n",
       "\n",
       "    var fmt_picker_span = $('<span/>');\n",
       "\n",
       "    var fmt_picker = $('<select/>');\n",
       "    fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content');\n",
       "    fmt_picker_span.append(fmt_picker);\n",
       "    nav_element.append(fmt_picker_span);\n",
       "    this.format_dropdown = fmt_picker[0];\n",
       "\n",
       "    for (var ind in mpl.extensions) {\n",
       "        var fmt = mpl.extensions[ind];\n",
       "        var option = $(\n",
       "            '<option/>', {selected: fmt === mpl.default_extension}).html(fmt);\n",
       "        fmt_picker.append(option)\n",
       "    }\n",
       "\n",
       "    // Add hover states to the ui-buttons\n",
       "    $( \".ui-button\" ).hover(\n",
       "        function() { $(this).addClass(\"ui-state-hover\");},\n",
       "        function() { $(this).removeClass(\"ui-state-hover\");}\n",
       "    );\n",
       "\n",
       "    var status_bar = $('<span class=\"mpl-message\"/>');\n",
       "    nav_element.append(status_bar);\n",
       "    this.message = status_bar[0];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.request_resize = function(x_pixels, y_pixels) {\n",
       "    // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n",
       "    // which will in turn request a refresh of the image.\n",
       "    this.send_message('resize', {'width': x_pixels, 'height': y_pixels});\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.send_message = function(type, properties) {\n",
       "    properties['type'] = type;\n",
       "    properties['figure_id'] = this.id;\n",
       "    this.ws.send(JSON.stringify(properties));\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.send_draw_message = function() {\n",
       "    if (!this.waiting) {\n",
       "        this.waiting = true;\n",
       "        this.ws.send(JSON.stringify({type: \"draw\", figure_id: this.id}));\n",
       "    }\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
       "    var format_dropdown = fig.format_dropdown;\n",
       "    var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n",
       "    fig.ondownload(fig, format);\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype.handle_resize = function(fig, msg) {\n",
       "    var size = msg['size'];\n",
       "    if (size[0] != fig.canvas.width || size[1] != fig.canvas.height) {\n",
       "        fig._resize_canvas(size[0], size[1]);\n",
       "        fig.send_message(\"refresh\", {});\n",
       "    };\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_rubberband = function(fig, msg) {\n",
       "    var x0 = msg['x0'] / mpl.ratio;\n",
       "    var y0 = (fig.canvas.height - msg['y0']) / mpl.ratio;\n",
       "    var x1 = msg['x1'] / mpl.ratio;\n",
       "    var y1 = (fig.canvas.height - msg['y1']) / mpl.ratio;\n",
       "    x0 = Math.floor(x0) + 0.5;\n",
       "    y0 = Math.floor(y0) + 0.5;\n",
       "    x1 = Math.floor(x1) + 0.5;\n",
       "    y1 = Math.floor(y1) + 0.5;\n",
       "    var min_x = Math.min(x0, x1);\n",
       "    var min_y = Math.min(y0, y1);\n",
       "    var width = Math.abs(x1 - x0);\n",
       "    var height = Math.abs(y1 - y0);\n",
       "\n",
       "    fig.rubberband_context.clearRect(\n",
       "        0, 0, fig.canvas.width, fig.canvas.height);\n",
       "\n",
       "    fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_figure_label = function(fig, msg) {\n",
       "    // Updates the figure title.\n",
       "    fig.header.textContent = msg['label'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_cursor = function(fig, msg) {\n",
       "    var cursor = msg['cursor'];\n",
       "    switch(cursor)\n",
       "    {\n",
       "    case 0:\n",
       "        cursor = 'pointer';\n",
       "        break;\n",
       "    case 1:\n",
       "        cursor = 'default';\n",
       "        break;\n",
       "    case 2:\n",
       "        cursor = 'crosshair';\n",
       "        break;\n",
       "    case 3:\n",
       "        cursor = 'move';\n",
       "        break;\n",
       "    }\n",
       "    fig.rubberband_canvas.style.cursor = cursor;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_message = function(fig, msg) {\n",
       "    fig.message.textContent = msg['message'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_draw = function(fig, msg) {\n",
       "    // Request the server to send over a new figure.\n",
       "    fig.send_draw_message();\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_image_mode = function(fig, msg) {\n",
       "    fig.image_mode = msg['mode'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.updated_canvas_event = function() {\n",
       "    // Called whenever the canvas gets updated.\n",
       "    this.send_message(\"ack\", {});\n",
       "}\n",
       "\n",
       "// A function to construct a web socket function for onmessage handling.\n",
       "// Called in the figure constructor.\n",
       "mpl.figure.prototype._make_on_message_function = function(fig) {\n",
       "    return function socket_on_message(evt) {\n",
       "        if (evt.data instanceof Blob) {\n",
       "            /* FIXME: We get \"Resource interpreted as Image but\n",
       "             * transferred with MIME type text/plain:\" errors on\n",
       "             * Chrome.  But how to set the MIME type?  It doesn't seem\n",
       "             * to be part of the websocket stream */\n",
       "            evt.data.type = \"image/png\";\n",
       "\n",
       "            /* Free the memory for the previous frames */\n",
       "            if (fig.imageObj.src) {\n",
       "                (window.URL || window.webkitURL).revokeObjectURL(\n",
       "                    fig.imageObj.src);\n",
       "            }\n",
       "\n",
       "            fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n",
       "                evt.data);\n",
       "            fig.updated_canvas_event();\n",
       "            fig.waiting = false;\n",
       "            return;\n",
       "        }\n",
       "        else if (typeof evt.data === 'string' && evt.data.slice(0, 21) == \"data:image/png;base64\") {\n",
       "            fig.imageObj.src = evt.data;\n",
       "            fig.updated_canvas_event();\n",
       "            fig.waiting = false;\n",
       "            return;\n",
       "        }\n",
       "\n",
       "        var msg = JSON.parse(evt.data);\n",
       "        var msg_type = msg['type'];\n",
       "\n",
       "        // Call the  \"handle_{type}\" callback, which takes\n",
       "        // the figure and JSON message as its only arguments.\n",
       "        try {\n",
       "            var callback = fig[\"handle_\" + msg_type];\n",
       "        } catch (e) {\n",
       "            console.log(\"No handler for the '\" + msg_type + \"' message type: \", msg);\n",
       "            return;\n",
       "        }\n",
       "\n",
       "        if (callback) {\n",
       "            try {\n",
       "                // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n",
       "                callback(fig, msg);\n",
       "            } catch (e) {\n",
       "                console.log(\"Exception inside the 'handler_\" + msg_type + \"' callback:\", e, e.stack, msg);\n",
       "            }\n",
       "        }\n",
       "    };\n",
       "}\n",
       "\n",
       "// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n",
       "mpl.findpos = function(e) {\n",
       "    //this section is from http://www.quirksmode.org/js/events_properties.html\n",
       "    var targ;\n",
       "    if (!e)\n",
       "        e = window.event;\n",
       "    if (e.target)\n",
       "        targ = e.target;\n",
       "    else if (e.srcElement)\n",
       "        targ = e.srcElement;\n",
       "    if (targ.nodeType == 3) // defeat Safari bug\n",
       "        targ = targ.parentNode;\n",
       "\n",
       "    // jQuery normalizes the pageX and pageY\n",
       "    // pageX,Y are the mouse positions relative to the document\n",
       "    // offset() returns the position of the element relative to the document\n",
       "    var x = e.pageX - $(targ).offset().left;\n",
       "    var y = e.pageY - $(targ).offset().top;\n",
       "\n",
       "    return {\"x\": x, \"y\": y};\n",
       "};\n",
       "\n",
       "/*\n",
       " * return a copy of an object with only non-object keys\n",
       " * we need this to avoid circular references\n",
       " * http://stackoverflow.com/a/24161582/3208463\n",
       " */\n",
       "function simpleKeys (original) {\n",
       "  return Object.keys(original).reduce(function (obj, key) {\n",
       "    if (typeof original[key] !== 'object')\n",
       "        obj[key] = original[key]\n",
       "    return obj;\n",
       "  }, {});\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.mouse_event = function(event, name) {\n",
       "    var canvas_pos = mpl.findpos(event)\n",
       "\n",
       "    if (name === 'button_press')\n",
       "    {\n",
       "        this.canvas.focus();\n",
       "        this.canvas_div.focus();\n",
       "    }\n",
       "\n",
       "    var x = canvas_pos.x * mpl.ratio;\n",
       "    var y = canvas_pos.y * mpl.ratio;\n",
       "\n",
       "    this.send_message(name, {x: x, y: y, button: event.button,\n",
       "                             step: event.step,\n",
       "                             guiEvent: simpleKeys(event)});\n",
       "\n",
       "    /* This prevents the web browser from automatically changing to\n",
       "     * the text insertion cursor when the button is pressed.  We want\n",
       "     * to control all of the cursor setting manually through the\n",
       "     * 'cursor' event from matplotlib */\n",
       "    event.preventDefault();\n",
       "    return false;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
       "    // Handle any extra behaviour associated with a key event\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.key_event = function(event, name) {\n",
       "\n",
       "    // Prevent repeat events\n",
       "    if (name == 'key_press')\n",
       "    {\n",
       "        if (event.which === this._key)\n",
       "            return;\n",
       "        else\n",
       "            this._key = event.which;\n",
       "    }\n",
       "    if (name == 'key_release')\n",
       "        this._key = null;\n",
       "\n",
       "    var value = '';\n",
       "    if (event.ctrlKey && event.which != 17)\n",
       "        value += \"ctrl+\";\n",
       "    if (event.altKey && event.which != 18)\n",
       "        value += \"alt+\";\n",
       "    if (event.shiftKey && event.which != 16)\n",
       "        value += \"shift+\";\n",
       "\n",
       "    value += 'k';\n",
       "    value += event.which.toString();\n",
       "\n",
       "    this._key_event_extra(event, name);\n",
       "\n",
       "    this.send_message(name, {key: value,\n",
       "                             guiEvent: simpleKeys(event)});\n",
       "    return false;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.toolbar_button_onclick = function(name) {\n",
       "    if (name == 'download') {\n",
       "        this.handle_save(this, null);\n",
       "    } else {\n",
       "        this.send_message(\"toolbar_button\", {name: name});\n",
       "    }\n",
       "};\n",
       "\n",
       "mpl.figure.prototype.toolbar_button_onmouseover = function(tooltip) {\n",
       "    this.message.textContent = tooltip;\n",
       "};\n",
       "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to  previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Pan axes with left mouse, zoom with right\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n",
       "\n",
       "mpl.extensions = [\"eps\", \"jpeg\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\"];\n",
       "\n",
       "mpl.default_extension = \"png\";var comm_websocket_adapter = function(comm) {\n",
       "    // Create a \"websocket\"-like object which calls the given IPython comm\n",
       "    // object with the appropriate methods. Currently this is a non binary\n",
       "    // socket, so there is still some room for performance tuning.\n",
       "    var ws = {};\n",
       "\n",
       "    ws.close = function() {\n",
       "        comm.close()\n",
       "    };\n",
       "    ws.send = function(m) {\n",
       "        //console.log('sending', m);\n",
       "        comm.send(m);\n",
       "    };\n",
       "    // Register the callback with on_msg.\n",
       "    comm.on_msg(function(msg) {\n",
       "        //console.log('receiving', msg['content']['data'], msg);\n",
       "        // Pass the mpl event to the overriden (by mpl) onmessage function.\n",
       "        ws.onmessage(msg['content']['data'])\n",
       "    });\n",
       "    return ws;\n",
       "}\n",
       "\n",
       "mpl.mpl_figure_comm = function(comm, msg) {\n",
       "    // This is the function which gets called when the mpl process\n",
       "    // starts-up an IPython Comm through the \"matplotlib\" channel.\n",
       "\n",
       "    var id = msg.content.data.id;\n",
       "    // Get hold of the div created by the display call when the Comm\n",
       "    // socket was opened in Python.\n",
       "    var element = $(\"#\" + id);\n",
       "    var ws_proxy = comm_websocket_adapter(comm)\n",
       "\n",
       "    function ondownload(figure, format) {\n",
       "        window.open(figure.imageObj.src);\n",
       "    }\n",
       "\n",
       "    var fig = new mpl.figure(id, ws_proxy,\n",
       "                           ondownload,\n",
       "                           element.get(0));\n",
       "\n",
       "    // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n",
       "    // web socket which is closed, not our websocket->open comm proxy.\n",
       "    ws_proxy.onopen();\n",
       "\n",
       "    fig.parent_element = element.get(0);\n",
       "    fig.cell_info = mpl.find_output_cell(\"<div id='\" + id + \"'></div>\");\n",
       "    if (!fig.cell_info) {\n",
       "        console.error(\"Failed to find cell for figure\", id, fig);\n",
       "        return;\n",
       "    }\n",
       "\n",
       "    var output_index = fig.cell_info[2]\n",
       "    var cell = fig.cell_info[0];\n",
       "\n",
       "};\n",
       "\n",
       "mpl.figure.prototype.handle_close = function(fig, msg) {\n",
       "    var width = fig.canvas.width/mpl.ratio\n",
       "    fig.root.unbind('remove')\n",
       "\n",
       "    // Update the output cell to use the data from the current canvas.\n",
       "    fig.push_to_output();\n",
       "    var dataURL = fig.canvas.toDataURL();\n",
       "    // Re-enable the keyboard manager in IPython - without this line, in FF,\n",
       "    // the notebook keyboard shortcuts fail.\n",
       "    IPython.keyboard_manager.enable()\n",
       "    $(fig.parent_element).html('<img src=\"' + dataURL + '\" width=\"' + width + '\">');\n",
       "    fig.close_ws(fig, msg);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.close_ws = function(fig, msg){\n",
       "    fig.send_message('closing', msg);\n",
       "    // fig.ws.close()\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.push_to_output = function(remove_interactive) {\n",
       "    // Turn the data on the canvas into data in the output cell.\n",
       "    var width = this.canvas.width/mpl.ratio\n",
       "    var dataURL = this.canvas.toDataURL();\n",
       "    this.cell_info[1]['text/html'] = '<img src=\"' + dataURL + '\" width=\"' + width + '\">';\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.updated_canvas_event = function() {\n",
       "    // Tell IPython that the notebook contents must change.\n",
       "    IPython.notebook.set_dirty(true);\n",
       "    this.send_message(\"ack\", {});\n",
       "    var fig = this;\n",
       "    // Wait a second, then push the new image to the DOM so\n",
       "    // that it is saved nicely (might be nice to debounce this).\n",
       "    setTimeout(function () { fig.push_to_output() }, 1000);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_toolbar = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var nav_element = $('<div/>')\n",
       "    nav_element.attr('style', 'width: 100%');\n",
       "    this.root.append(nav_element);\n",
       "\n",
       "    // Define a callback function for later on.\n",
       "    function toolbar_event(event) {\n",
       "        return fig.toolbar_button_onclick(event['data']);\n",
       "    }\n",
       "    function toolbar_mouse_event(event) {\n",
       "        return fig.toolbar_button_onmouseover(event['data']);\n",
       "    }\n",
       "\n",
       "    for(var toolbar_ind in mpl.toolbar_items){\n",
       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
       "\n",
       "        if (!name) { continue; };\n",
       "\n",
       "        var button = $('<button class=\"btn btn-default\" href=\"#\" title=\"' + name + '\"><i class=\"fa ' + image + ' fa-lg\"></i></button>');\n",
       "        button.click(method_name, toolbar_event);\n",
       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
       "        nav_element.append(button);\n",
       "    }\n",
       "\n",
       "    // Add the status bar.\n",
       "    var status_bar = $('<span class=\"mpl-message\" style=\"text-align:right; float: right;\"/>');\n",
       "    nav_element.append(status_bar);\n",
       "    this.message = status_bar[0];\n",
       "\n",
       "    // Add the close button to the window.\n",
       "    var buttongrp = $('<div class=\"btn-group inline pull-right\"></div>');\n",
       "    var button = $('<button class=\"btn btn-mini btn-primary\" href=\"#\" title=\"Stop Interaction\"><i class=\"fa fa-power-off icon-remove icon-large\"></i></button>');\n",
       "    button.click(function (evt) { fig.handle_close(fig, {}); } );\n",
       "    button.mouseover('Stop Interaction', toolbar_mouse_event);\n",
       "    buttongrp.append(button);\n",
       "    var titlebar = this.root.find($('.ui-dialog-titlebar'));\n",
       "    titlebar.prepend(buttongrp);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._root_extra_style = function(el){\n",
       "    var fig = this\n",
       "    el.on(\"remove\", function(){\n",
       "\tfig.close_ws(fig, {});\n",
       "    });\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._canvas_extra_style = function(el){\n",
       "    // this is important to make the div 'focusable\n",
       "    el.attr('tabindex', 0)\n",
       "    // reach out to IPython and tell the keyboard manager to turn it's self\n",
       "    // off when our div gets focus\n",
       "\n",
       "    // location in version 3\n",
       "    if (IPython.notebook.keyboard_manager) {\n",
       "        IPython.notebook.keyboard_manager.register_events(el);\n",
       "    }\n",
       "    else {\n",
       "        // location in version 2\n",
       "        IPython.keyboard_manager.register_events(el);\n",
       "    }\n",
       "\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
       "    var manager = IPython.notebook.keyboard_manager;\n",
       "    if (!manager)\n",
       "        manager = IPython.keyboard_manager;\n",
       "\n",
       "    // Check for shift+enter\n",
       "    if (event.shiftKey && event.which == 13) {\n",
       "        this.canvas_div.blur();\n",
       "        event.shiftKey = false;\n",
       "        // Send a \"J\" for go to next cell\n",
       "        event.which = 74;\n",
       "        event.keyCode = 74;\n",
       "        manager.command_mode();\n",
       "        manager.handle_keydown(event);\n",
       "    }\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
       "    fig.ondownload(fig, null);\n",
       "}\n",
       "\n",
       "\n",
       "mpl.find_output_cell = function(html_output) {\n",
       "    // Return the cell and output element which can be found *uniquely* in the notebook.\n",
       "    // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n",
       "    // IPython event is triggered only after the cells have been serialised, which for\n",
       "    // our purposes (turning an active figure into a static one), is too late.\n",
       "    var cells = IPython.notebook.get_cells();\n",
       "    var ncells = cells.length;\n",
       "    for (var i=0; i<ncells; i++) {\n",
       "        var cell = cells[i];\n",
       "        if (cell.cell_type === 'code'){\n",
       "            for (var j=0; j<cell.output_area.outputs.length; j++) {\n",
       "                var data = cell.output_area.outputs[j];\n",
       "                if (data.data) {\n",
       "                    // IPython >= 3 moved mimebundle to data attribute of output\n",
       "                    data = data.data;\n",
       "                }\n",
       "                if (data['text/html'] == html_output) {\n",
       "                    return [cell, data, j];\n",
       "                }\n",
       "            }\n",
       "        }\n",
       "    }\n",
       "}\n",
       "\n",
       "// Register the function which deals with the matplotlib target/channel.\n",
       "// The kernel may be null if the page has been refreshed.\n",
       "if (IPython.notebook.kernel != null) {\n",
       "    IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n",
       "}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAYAAAB91L6VAAAgAElEQVR4Xu2dC9ytU534v8dd6JU6Mh25lgozhUKjPwYzNYlOUoNuuk2RxiSV9E+lSFLTRQ1lItTkVtLMFErl1j+XUhyaSCqXOm4HRRw6/8/vWFu7bd/3fvfzrLW/z+fzfg7vc1m/9f2t/X73Ws961jMHNwlIQAISkIAEJk5gzsRLtEAJSEACEpCABFDANgIJSEACEpBABQQUcAXQLVICEpCABCSggG0DEpCABCQggQoIKOAKoFukBCQgAQlIQAHbBiQgAQlIQAIVEFDAFUC3SAlIQAISkIACtg1IQAISkIAEKiCggCuAbpESkIAEJCABBWwbkIAEJCABCVRAQAFXAN0iJSABCUhAAgrYNiABCUhAAhKogIACrgC6RUpAAhKQgAQUsG1AAhKQgAQkUAEBBVwBdIuUgAQkIAEJKGDbgAQkIAEJSKACAgq4AugWKQEJSEACElDAtgEJSEACEpBABQQUcAXQLVICEpCABCSggG0DEpCABCQggQoIKOAKoFukBCQgAQlIQAHbBiQgAQlIQAIVEFDAFUC3SAlIQAISkIACtg1IQAISkIAEKiCggCuAbpESkIAEJCABBWwbkIAEJCABCVRAQAFXAN0iJSABCUhAAgrYNiABCUhAAhKogIACrgC6RUpAAhKQgAQUsG1AAhKQgAQkUAEBBVwBdIuUgAQkIAEJKGDbgAQkIAEJSKACAgq4AugWKQEJSEACElDAtgEJSEACEpBABQQUcAXQLVICEpCABCSggG0DEpCABCQggQoIKOAKoFukBCQgAQlIQAHbBiQgAQlIQAIVEFDAFUC3SAlIQAISkIACtg1IQAISkIAEKiCggCuAbpESyIzA9sB3geWBBzOL3XAlUFsCCri2qTGwTAh8D/hb4IGWeF8PnJJJHSLMG4C1WgT7NeBVgALOKJGGmg8BBZxProy0ngRCwBcC/7ee4fUdVQj4Q8Bxbc5QwH1j9EAJ9E9AAffPyiMl0I5ALwGH0DYFtk295KcClwJvBE5LF3wN8DZgfeCPwFeA/dO+pwNHAc9JvdOvA+8A/gDE5/cDwOuAGeAu4ETgYGAF4BPAbsAqwO3Ax4BPd0jjIAI+AVgOeGXTtZo5hLD/O/Wco64R538BfwJ2BZbYlCQggYc/GG4SkMDwBHoJeKXUQ/5hEmf8ey5wQCryDcCHgT2BuNbKwObA94EnANcAhwOfBVZLcv5lEvjfA8enIfBfA2sATwN+kPbvC8Qxt6Xh5ScBP5qAgKOI+JKwX6rLW4D4krEFsGh41J4pgbIITKuAP5W+ia8LbAZcMaa0Rk8kejLPAA5MPZDGpbvtG1PxXqYCAiHNrVPPtbn46LFem36xHnAZEJKMnuvfNd1rvQqIHmX0clu3kPTLgOc27dgGOA94DPA84Ix0nzbiuK/puBBeDIvHvegQ8uIebKIHvGbLvewXpy8CrUPQvXrAjaJOB56SfiLWcX3OKkizRUpg/ASmVcAxHHh96pnMH+MfhmemP2DvTj2NGAJsbN32jT+zXnFSBHr1gBtxfC71Sp8PnNMUXAh5D+AbbQL+9yTQe5v2xWd2RWBD4KY0/ByybXyR/GDqYS+bvgyGwGMY++Ik5Ms7gBnnEHSjiJicdhHwn8Bek0qI5UggFwLTKuBGfuKPTrOAo9fyEeCxQPwBi6G/xn26QXIaPYT4tt8s4Mb53fYNUobH1oNAPwLeOUno1HRfNIZi435tbAuAL6T7s601eh+wHbBDH1UNKcdQb7TZucA9TeesChwKvBxYewwCjvvIcZ2XNF3r50DUrzEZLe5JR68/htKj3Fd0+JLRR9U8RAJlElDAfxbw6ulZxxcCt6T7b3G/LIb/oqcxyKaAB6GV97G9BLxBEtE+6ctcTE66P8krJiP9c5p9HL3gkFUMLTfuAa8D/AR4TxqmjiHmEF8I/ExgSyDuMcdEp5i8FfeT4/bK49Ow+N3AT9Pkp0PSfdi47dJuG6QHHD3uT6YyfpHE/9H05TUEHH9X4hGmmKi1S+rhxz3siDtGntwkIIH0QZlmEM094BBvDJXFBJfGFpNa9k733G5Mf+za8YrJMs2bAp6eVtXpOeDocUZPMYZ+45iYGxBbtKn4YncMcET6XdynjVnPca84hpujHTaOj+HjmKQV95lDzvFlMPbHUHP0jI9ME69ihnH0Qt8LfCtJL2QY14z7vyHimJcQsh5VwCHWo1PPNp5/jrrEfeLG41gHpS8WIdw7U2HBIu4Dx7B0873q6Wkp1lQCLQTsAf+5BxzDhNHTiD8Qo24KeFSCni8BCUigcAIK+M8CfhxwdZpR+u2U92el37WuctSrWSjgXoTcLwEJSGDKCUyrgI8FoscbS+/FAgUxYSUel4h7b/E4SNxDi3Vv47GRmKQV99f62WK4OlYTCpnHsN/v0z2wH6eh7E77+rm2x0hAAhKQQEEEplXABaXQqkhAAhKQQI4EFHCOWTNmCUhAAhLInoACzj6FVkACEpCABHIkoIBzzJoxS0ACEpBA9gSmTsArrLDCkrlzY6EgNwlIQAISyJnATTfdFE+oxCpwWW5TJ+B58+YtufHGWFPDTQISkIAEciYwZ86cWJim0/Kqta+aAq59igxQAhKQgATaEVDAmbULe8CZJcxwJSABCXQgoIAzaxoKOLOEGa4EJCABBVxGG1DAZeTRWkhAAhKwB5xZG1DAmSXMcCUgAQnYAy6jDSjgMvJoLSQgAQnYA86sDSjgzBJmuBKQgATsAZfRBhRwGXm0FhKQgATsAWfWBhRwZgkzXAlIQAL2gMtoAwq4jDxaCwlIQAL2gDNrAwo4s4QZrgQkIAF7wGW0AQVcRh6thQQkIAF7wJm1AQWcWcIMVwISkIA94DLagAIuI4/WQgISkIA94MzawKgCvmSTS5bWeMsFW2ZWc8OVgAQkUBYBBZxZPhVwZgkzXAlIQAIOQZfRBhRwGXm0FhKQgATsAWfWBhRwZgkzXAlIQAL2gMtoAwq4jDxaCwlIQAL2gDNrAwo4s4QZrgQkIAF7wGW0AQVcRh6thQQkIAF7wJm1AQWcWcIMVwISkIA94DLagAIuI4/WQgISkIA94MzagALOLGGGKwEJSMAecBltQAGXkUdrIQEJSMAecGZtQAFnljDDlYAEJGAPuIw2oIDLyKO1kIAEJGAPOLM2oIAzS5jhSkACErAHXEYbUMBl5NFaSEACErAHPHgbWAn4CrAxcB+wENgHuK7lUusBvwCubPr9S9Pv4lcvAo4Clk3H7A3c3SscBdyLkPslIAEJ5EFAAQ+epxDwDsA3gSXAfsDuwPZtBHwFsHqbIlZNIt4O+BlwdJL5O3qFo4B7EXK/BCQggTwIKODR8/Rs4HQgerzNW/x/JwG/DHg98IJ0QvSmzwHW7hXOqAK++MkX88AtD7DJ6Zswd/7cXsW5XwISkIAEZomAAh4d7EnAHcD+bQR8bZJwDDOfCRwGPAS8HdgIeFM65zHAPcCKwIPdQhpFwLeeeSsLXrJg6eXnLD+HjU/dWAmPnn+vIAEJSGAoAgp4KGyPnHQwsAuwI3Bvy6VCpjPpHvEawCnAucCRAwr4ACB+lm4zMzPzFi1aNHDUId+rX341SxbHqPnDmxIeGKMnSEACEhgbAQU8PMoDgT2AnYB+jLgnsFcS9kSHoNvJVwkPn3jPlIAEJDAOAgp4OIrRI31Fku+dHS6xJhD7Fqeh5ZOBa4BDgNXSJKxtmyZh/REIqXfdBh2C7iZfJdyLtvslIAEJzB4BBTw425go9Rvg+nTfNq5wP7AVcChwM3AMsFv6/7jnuxxwXhJsHBvbrmk4OvZdBbwGuKtXOIMK+KI1L2LxrfEdoPu2/Nzl2WbhNr0Oc78EJCABCYyJgAIeE8hJXWZQAdsDnlRmLEcCEpDAYAQU8GC8Kj96UAFHwN4DrjxtBiABCUjgUQQUcGaNYhgBd5Kws6AzS77hSkACRRFQwJmlc1gBNyTsc8CZJdxwJSCBYgko4MxSO4qAo6quhJVZwg1XAhIoloACziy1owr4kk0uWVrjLRdsmVnNDVcCEpBAWQQUcGb5VMCZJcxwJSABCXQgoIAzaxoKOLOEGa4EJCABBVxGG1DAZeTRWkhAAhKwB5xZG1DAmSXMcLMj4DyJ7FKWbcAKOLPUKeDMEma42RHwSYHsUpZtwAo4s9Qp4MwSZrhZEfCd2VmlK/tgFXBmKVTAmSXMcLMh4Duzs0lVMYEq4MxSqYAzS5jhZkHA9dKzSFNxQSrgzFKqgDNLmOHWnoBvDKt9iooNUAFnlloFnFnCDLf2BHxndu1TVGyACjiz1CrgzBJmuLUnYA+49ikqNkAFnFlqFXBmCTPcLAh4DziLNBUXpALOLKUKOLOEGW42BJwFnU2qiglUAWeWSgWcWcIMNysCPgecVbqyD1YBZ5ZCBZxZwgw3OwKuhJVdyrINWAFnlrpRBZxZdQ1XAhKQQLEEFHBmqVXAmSXMcCUgAQl0IKCAM2saCjizhBmuBCQgAQVcRhtQwGXk0VpIQAISsAecWRtQwJklzHAlIAEJ2AMuow0o4DLyaC0kIAEJ2APOrA0o4MwSZrgSkIAE7AGX0QYUcBl5tBYSkIAE7AFn1gYUcGYJM1wJSEAC9oDLaAMKuIw8WgsJSEAC9oAzawMKOLOEGa4EJCABe8BltAEFXEYerYUEJCABe8CZtQEFnFnCDFcCEpCAPeAy2oACLiOP1kICEpCAPeDB28BKwFeAjYH7gIXAPsB1LZf6a+AzwJrAg8AlwFvSOXHoEuAq4KF03luBC3qFo4B7EXK/BCQggTwIKODB8xQC3gH4ZpLofsDuwPYtl3oqsDLwU2BZ4MvANcD703Eh4McBiwYJQQEPQstjJSABCdSXgAIePTfPBk4H1utxqQOBTYG9FfDo0L2CBCQggdwJKODRM3gScAewf5dLrQJcDrwb+FqTgH8ELAd8B3gv8Ide4dgD7kXI/RKQgATyIKCAR8vTwcAuwI7AvR0utQLwVeB64F+ajlkH+DUQcj4GuAfYt801DgDiZ+k2MzMzb9GigUatR6uhZ0tAAhKQwKwQUMDDY40h5T2Anbrcx10eOBW4DfjndM+4XYnPBT4HxMStrps94F6E3C8BCUggDwIKeLg8RY/0FUm+d3a4RAwtn5Lk/IYW+cbkq/tTr3kZ4OPAGsCre4WjgHsRcr8EJCCBPAiUIODzRkQds5FjCLnfbW3gN2lIOYaNYwuZbgUcCtychpRD0CenWdBRRmwXpUeRosd7bJJyiDruBcc95LiXbA+4FwT3S0ACEiiAQAkC/tOIeQg5xmNCWWz2gLNIk0FKQAIS6EmgBAFv17OWnQ84CPgHBTwCQU+VgAQkIIGhCJQg4KEqnk46Pt13tQc8CkXPlYAEJCCBgQlMu4A/DMwHnjEwuYpOcAi6IvAWKwEJSGDMBKZdwGPGOfuXU8Czz9gSJCABCUyCgAKeBOUxlqGAxwjTS0lAAhKokIACrhD+MEUr4GGoeY4EJCCB+hFQwPXLSdeIFHBmCTNcCUhAAh0IlCjg1YDNgPNLzLoCLjGr1kkCEphGAiUION65e19T8mJFqotzerZ3kIangAeh5bESkIAE6kugBAH/GHgB8LuEWQHXt70ZmQQkIAEJJAIlCPjnqbf7fOC6tCazPWCbuAQkIAEJ1JpACQJ+PPAtIN6v+49AvAJQAde62RmcBCQgAQmUIODI4qrAWcAWwEeAD3oP2MYtAQlIQAJ1JlCKgIPxCun9uy9Or/nLZn3nQRqIk7AGoeWxEpCABOpLoCQBB+V4uf1ngbgfvH59sQ8fmQIenp1nSkACEqgTgdIEXCe2sxKLAp4VrF5UAhKQwMQJKOCJIx+tQAU8Gj/PloAEJFAXAiUIeDfgbOAPdYE6m3Eo4Nmk67UlIAEJTI5ACQJ+KL3PN54HLn5TwMWn2ApKQAJTQqAEAf8JeDqggKek0VpNCUhAAiUQUMCZZdEecGYJM1wJSEACHQgo4MyahgLOLGGGKwEJSEABl9EGFHAZebQWEpCABOwBZ9YGFHBmCTNcCUhAAvaAy2gDCriMPFoLCUhAAvaAM2sDCjizhBmuBCQgAXvAZbQBBVxGHq2FBCQggRJ6wK8BvgrcMw3pVMDTkGXrKAEJTAOBEgQ8DXl6pI4KeKrSbWUlIIGCCSjgzJKrgDNLmOFKQAISKPgecGvVVgT+Cri1xBc0KGA/yxKQgATKIFBSD3h14PPAfGCZlJ5fApcDl6V/47/vyjl1Cjjn7Bm7BCQggT8TKEnAxwGvA5YAvwUeC6ySqhq/a2zXA5cCe+XYEBRwjlkzZglIQAKPJlCSgG8GZoDtUo83arsRsAWwefp3s3RMCHnZHBuEAs4xa8YsAQlIoGwB3w18B3hJj0RvmGR86pANYiXgK8DGwH3AQmAf4Lo213sRcFSS/ZXA3kDEGVu3fR1DU8BDZs3TJCABCdSMQEk94B+koedeAh41BSHgHYBvpuHu/YDdge1bLrwq8IvUI/8ZcHQS9juAbvu6xqeAR02f50tAAhKoB4GSBPwm4MPAOsDvJ4j32cDpwHotZb4MeD3wgvT76DGfA6wNdNungCeYPIuSgAQkUBWBkgQ8B/h+GhKOCVYPTAjqScAdwP4t5b093YOOLwaxPSat1hWPScWxcX+63b4Hu8VtD3hCWbUYCUhAArNMoCQBB6qYhHVJEt0H0zDxbIr4YGAXYEfg3lkS8AFA/CzdZmZm5i1atGiWm4WXl4AEJCCB2SZQkoBjCDjuy0bPMnrDMdP5j8BP0mNH8SxwPH4U92PHsR0I7AHsBLQzokPQ46DsNSQgAQkUSqAkAX8NeHGa6HQVsBrw1KbHjRrPAsf94R8BfzdCTqNH+ook3zs7XCfKj0lY2ybpxySs+EIQ4u62r2tYDkGPkDVPlYAEJFAjAiUJOO7DhmRjUlSsgBVb3HeNZ3/jWeD4ffw0esjDPgcck6h+A8SCHo03MN0PbAUcCsTzyMek8ncFjgSWA+JLQby5qbESV7d9HZuIAq7Rp8dQJCABCYxAoCQBxzDwuWmGcTck8QhQSPmCEbhVdqoCrgy9BUtAAhIYK4GSBPzdNBFq57ESqtnFFHDNEmI4EpCABIYkUJKAXw58AYiVrn43JI/an6aAa58iA5SABCTQF4GSBByPIH0aWB+I+6udJkf1BaauByngumbGuCQgAQkMRqAkAT8E3ATMA24EPg58NU2YGoxKjY9WwDVOjqFJQAISGIBASQKOSVjxCsLG1njsKF5NGM8AN54Djn9vG4BRrQ5VwLVKh8FIQAISGJpASQIOCI03HTVePxj/Pi7RaX4n8K/TUPXQ4Ko6UQFXRd5yJSABCYyXQGkCbkcnVsiK54Ab7wWOZ4FDysM+BzzeDAx4NQU8IDAPl4AEJFBTAtMg4Hbo441J0QvOblPA2aXMgCUgAQm0JTCtAs62OSjgbFNn4BKQgAT+gkBJAn4dsDfwAeA7peZZAZeaWeslAQlMG4GSBBxvQooXLMRazZ1mOceazF3ft1v3BqCA654h45OABCTQH4GSBNx4QcJ2Xaoerwg8CngrcFZ/iOp1lAKuVz6MRgISkMCwBEoScLzq7zTgVV1grJSWqTwbiKUrs9sUcHYpM2AJSEACbQmUJOCFwA+BXXrk+kJgLeApObYJBZxj1oxZAhKQwKMJlCTgeBVhvGbwScADXZL9n8CLgNVybBAKOMesGbMEJCCBsgW8L3A08H+Bw7sk+9vA3wKPybFBKOAcs2bMEpCABMoW8MrAT9MSkyHhI9ok/InAdekFDRvn2CAUcI5ZM2YJSEACZQs4avc3QPRwH59kfAxwEXAPEML9EPAs4N+AA3NsEAo4x6wZswQkIIHyBRw13Aj4MhAvYmh+AUPsm5NeWRj3irN8I5IC9mMsAQlIoAwCJU3Cas5IiPafgHjuN0S8JhCvKzwHeG96X3CWGVTAWabNoCUgAQk8ikCpAi421Qq42NRaMQlIYMoIKODMEq6AM0uY4UpAAhLoQEABZ9Y0FHBmCTNcCUhAAgq4jDaggMvIo7WQgAQkYA84szaggDNLmOFKQAISsAdcRhtQwGXk0VpIQAISsAecWRtQwJklzHAlIAEJ2AMuow0o4DLyaC0kIAEJ2APOrA0o4MwSZrgSkIAEpqgHvA6wPXBiiVlXwCVm1TpJQALTSKDEHvBLgVOBZUtMqAIuMavWSQISmEYCCjizrCvgzBJmuBKQgAQKHoK+vqVujwHm2gO2zUtAAhKQQJ0JlNAD/vcWwBsAOyngOjc7Y5OABCQggRIE3JrFSdwD/hSwK7AuEO8WvqJNU3otsH/T79cGzgd2A9YDfgFc2bQ/4o7fdd0cgu5FyP0SkIAE8iBQqoBPA5aZxRRsC8TQ94XA/A4Cbi3+KuB9wBlJwCHt1QeNUQEPSszjJSABCdSTQIkCXhlYE/jVBJDf0KeAtwK+AcwDFivgCWTGIiQgAQnUnECJAp4k8n4F/DngHuDtKbgYgr429ZzjcakzgcOAh3oFbw+4FyH3S0ACEsiDgAIeLU/9CHgV4BZga+DqVNyKwAywEFgDOAU4FziyTTgHAPGzdJuZmZm3aNGi0aL2bAlIQAISqJxACQKOSU1nA3+ogGY/At4beBPw3C7x7QnsBezSqw72gHsRcr8EJCCBPAiUIOAYtn0G8PMKkPcj4AuALwLHNcUX96jvTPeDozd8MnANcEivOijgXoTcLwEJSCAPAiUI+E/A0ycs4GOBnYG1gNvT/d2nJMmeBcRPbE8DLgOelI5ptIrotR+a7vkuB5wHHAjc36vZKOBehNwvAQlIIA8CCjiPPD0SpQLOLGGGKwEJSKADAQWcWdNQwJklzHAlIAEJKOAy2oACLiOP1kICEpCAPeDM2oACzixhhisBCUjAHnAZbUABl5FHayEBCUjAHnBmbUABZ5Yww5WABCRgD7iMNqCAy8ijtZCABCRQQg/4NcBXW56zLTazCrjY1FoxCUhgygiUIOCpSpkCnqp0W1kJSKBgAgo4s+Qq4MwSZrgSkIAECr4H3KjaCsADpWdaAZeeYesnAQlMC4GSesAfAvYFXgzECxCat1izOd67e1PuiVXAuWfQ+CUgAQk8TKAkAV8EbAI8seWlBkcBb0sJvxR4FXBtrg1AAeeaOeOWgAQk8JcEShLwzcB1wLZNVdwU+CmwBLgDeDzwK+Cvgd/n2BgUcI5ZM2YJSEACjyZQkoDjVX7xOFK83L6xHQ4cBPwr8CngPcAHgYOBI3JsEAo4x6wZswQkIIGyBXwX8F1gflM1fwKsnYalHwSWST3guBe8dY4NQgHnmDVjloAEJFC2gC8H1gDWT9V8cpLtqcAeTVX/OrAdsHqODUIB55g1Y5aABCRQtoAbQ8tvBj4PHAO8EXgDcHxT1b8E7A6smGODUMA5Zs2YJSABCZQt4Oj9Xg3MBe4GHgv8AVgXuLOp6t8Hngo8KccGoYBzzJoxS0ACEihbwFG7eAwper9xf/de4E1A9HgbWyzWcXsS9VY5NggFnGPWjFkCEpBA+QJu1PAx6Vngh1qq/BLgDOCTTc8GZ9UuFHBW6TJYCUhAAh0JlPQYUj9p3hF4GfAV4Hv9nFC3YxRw3TJiPBKQgASGIzBtAh6OUo3OUsA1SoahSEACEhiBgAIeAV4VpyrgKqhbpgQkIIHxE1DA42c6q1dUwLOK14tLQAISmBgBBTwx1OMpSAGPh6NXkYAEJFA1AQVcdQYGLF8BDwjMwyUgAQnUlIACrmliOoWlgDNLmOFKQAIS6EBAAWfWNBRwZgkzXAlIQAIKuIw2oIDLyKO1kIAEJGAPOLM2oIAzS5jhSkACErAHXEYbUMBl5NFaSEACErAHnFkbUMCZJcxwJSABCUxBD/hpwP+WnmkFXHqGrZ8EJDAtBErqAf8pvQf4R8Bl6edS4JclJVMBl5RN6yIBCUwzgZIEvAB4CrBcS0LvbBJyiPliYOGISf8UsCuwLrAZcEWb620PfLOlV/5c4L507OuBg4BlgPOAfYHFveJSwL0IuV8CEpBAHgRKEnAQPyxJ7TbgemBlYCNgJWBJU0pCmMcA/wFEz3nQbdt0/QuB+V0E/AngWW0uvj5wEbA58Dvg68DZwGd6BaKAexFyvwQkIIE8CJQk4H2Ao4GDgY82iTV6xC8GPgQ07hNHTzl6njFc/aIkwWEydsOQAn4HsCHw5lToC1Pcz+sVhALuRcj9EpCABPIgUJKArwL+CDy7A/roBX8DeCzwT6mn/M/Aj4GtgAeHSFkvAUd51wIPAccDn01lfBq4Gfhw+v+NgW8B6/SKQQH3IuR+CUhAAnkQKEnAcW/1a8BeXdDPTZOy3plkGEPEbwXeCHxhiJR1E3CIfg5wF7A28D+pF34qMIiADwDiZ+k2MzMzb9GiRUOE6ikSkIAEJFAnAiUJ+Ebgt116wA3u/w2sAcSEqNWBXwP/D/iHIRLTTcCtl3s38KQkfIegh4DtKRKQgARKIlCSgE8AXpVE+p0uSToD2Ck6k+mY89O94ScOkdhuAv6rdG85JnmtloaYY9JX9LQ3AGICV/MkrHPSPeyuYTgEPUSWPEUCEpBADQmUJOCnp9nI9wKvBv6rDe9VgJ+n+8Ahxdi+DLwUWHGA/BwL7AysBdwO3JMegToOOCv97AfExLC4txwTwU4DPtA0GzuGveMxpNi+lyZk+RjSAEnwUAlIQAI5EyhJwJGHPYHoCYfw4tnak4FYjOP3SZDvBeIRotj39ylx8azu1sDjckikPeAcsmSMEpCABHoTKE3AUeNt0vO98fxv87O/sS8mRcVM6e2SmON3MYwcv4sedO03BVz7FBmgBCQggb4IlCTgWJHqyjTkG8/4xqNGu6VJWWumZSovAD6YjgtAWyQRf67pmdy+wFV1kAKuirzlSkACEhgvgZIEHJOdYjbz3w6IKCZfhbBvGfC8Sg5XwJVgt1AJSEACYydQkoDj4diYePXKsVOq0QUVcKNOjpUAAByCSURBVI2SYSgSkIAERiBQkoC/n1ac2mEEHrU/VQHXPkUGKAEJSKAvAiUJOHq+8YztpulRo74A5HaQAs4tY8YrAQlIoD2BkgQcC2vEQhfPAHZJbysqLu8KuLiUWiEJSGBKCZQk4HjhQbzgYF5aGONjQKy7/LOScquAS8qmdZGABKaZQEkCjklY8QKExtZ4BjhmN1+eHje6DIifeF9wlpsCzjJtBi0BCUjgUQRKEnBULt6xG8/2xhrLjX8bK1w1L8oRL2BYP8f2oIBzzJoxS0ACEng0gdIE3C7H6yUZN4Qc7wsOKS+bY4NQwDlmzZglIAEJTKeA2+U9XnwfveDsNgWcXcoMWAISkEBbAtPQAy4q9Qq4qHRaGQlIYIoJlCjg+cA/Ak8A7gTO7PBqwizTroCzTJtBS0ACEngUgZIEHPd0z0jPAMdbj2KLiVdfBF6X/j+kfAcQ60ZnuSngLNNm0BKQgASKFvA7gSOAXwKHpdcMfju9H7gh4Pj3k8CL0zuBs2sSCji7lBmwBCQggbYESuoBXwXE5Kp4r28syBFb9HRPaOoBx2pZvwWOAd6WY5tQwDlmzZglIAEJPJpASQK+D/heuv/bqGmrgOP38crC5dJ7grNrEwo4u5QZsAQkIIHie8CxEtZ5wG5NNW0n4Fiecpu0ZGV2zUIBZ5cyA5aABCRQvIB/ADwR2KCHgL8EvBRYKcc2oYBzzJoxS0ACEih7CPrdwIeAN6bXEkZt2/WAY5g6lqpsXjc6m7ahgLNJlYFKQAIS6EqgpHvAsbxkvPkoJlq9BjiljYDXBq4DFqTlKbNrHgo4u5QZsAQkIIHih6CjgtumRTdWAS5O93rjnu8bgE2BzwDPAj4IvD/HNqGAc8yaMUtAAhIoewi6UbsQ7IlJuPG75rcgxQId1wObpXcGZ9cmFHB2KTNgCUhAAlPRA25UchlgjzTZKmS7JnAX8D/Ae9OzwFk2CQWcZdoMWgISkMCjCJR0D3gq0quApyLNVlICEpgCAgo4syQr4MwSZrgSkIAEOhBQwJk1DQWcWcIMVwISkIACbktgZ+CZwOG5tBAFnEumjFMCEpBAdwLT3gM+Hng1EK8yzGJTwFmkySAlIAEJ9CSggBVwz0biARKQgAQkMH4CJQg4nusddnsCEIt22AMelqDnSUACEpDAUARKEHCs9zzKFgt1KOBRCHquBCQgAQkMTKAEAa87cK0ffcKvxnCNiVzCe8ATwWwhEpCABGadQAkCnnVIbQr4FLArEPKPlbauaHPMDsARwKppOcz/Bg5KL4hYD/gFcGXTefGKxPhd100B9yLkfglIQAJ5EFDAw+UpXvoQ954vBOZ3EHCIOZa/jOPi3cPfBo4DTgBCwCHt1QctXgEPSszjJSABCdSTgAIeLS83dBFw65WPBm5Lb2FSwKNx92wJSEAC2RNQwKOlsF8Br5V6vC8CLks94GvT72IC2JnAYcBDvcKxB9yLkPslIAEJ5EFAAY+Wp34E/FjgO8B/Ah9Pxa0IzAALgTWAU4BzgSPbhHMAED9Lt5mZmXmLFi0aLWrPloAEJCCBygko4NFS0EvAqwFnp9cgfqhLUXsCewG79ArHHnAvQu6XgAQkkAcBBTxanroJOGY/h3zj59CWYuL9xHcCi4HoDZ8MXAMc0iscBdyLkPslIAEJ5EFAAQ+Xp2OBeJFD3Nu9HbgHeEqa5XwWED/vSROuFjQVcVq617tbknLc810OOA84ELi/VzgKuBch90tAAhLIg4ACziNPj0SpgDNLmOFKQAIS6EBAAWfWNBRwZgkzXAlIQAIKuC2BJwN/nX4+kkMrUcA5ZMkYJSABCfQmMC094HjkpyHaxr+bAvGI0Jx073Xl3riqP0IBV58DI5CABCQwDgIlCziWi3wH8ExgXhOs3wE/BX6S/o3/jhnID44D6GxfQwHPNmGvLwEJSGAyBEoWcDwiFI8CfRWIVafeBDwOeCvw5cngHX8pCnj8TL2iBCQggSoIlCzgeE9wPCr0zQQ2XogQz+P+K/At4M3AzVVAH6VMBTwKPc+VgAQkUB8CJQt4a+B/04IXzcTjLUXxVqINgXcB8UxvNpsCziZVBioBCUigK4GSBdyt4suk9ZXfD1wKvKGfd/HWoS0p4DpkwRgkIAEJjE5gWgXcILcBcAywDbDK6Dhn/woKePYZW4IEJCCBSRAoVcDxaNGSAQC+CjhpgOMrO1QBV4begiUgAQmMlUCpAr4bOAI4fKy0anAxBVyDJBiCBCQggTEQKFXAMQP6BOB1Y2BUq0so4Fqlw2AkIAEJDE1AAQ+NrpoTFXA13C1VAhKQwLgJKOBxE53l6yngWQbs5SUgAQlMiIACnhDocRWjgMdF0utIQAISqJZACQJ+C3BV+rk94fQecLXtytIlIAEJSKAHgRIE/FBTHRcmEe+YlpuMBTayW26yW87sAfuZloAEJFAGgRIEfA6wObBGS0oazwHH248uSyteNf69Ldf0KeBcM2fcEpCABP6SQAkCbtRoXWCLpp+Q8hOaqtu8MMevgfVzbAwKOMesGbMEJCCBRxMoScDt8vvkFiGHoNdMq2Qtm2ODUMA5Zs2YJSABCUyfgNvlfF6S8lk5NggFnGPWjFkCEpCAAs6+DSjg7FNoBSQgAQksJVD6EHRxaVbAxaXUCklAAlNKQAFnlngFnFnCDFcCEpBABwLTIuDlgU2AzYBnAc8Ets+xVSjgHLNmzBKQgASm4x7wTJJsiLYh3KcDIeGlw+7OgvajIAEJSEACVRMooQe8a5NwQ7rxPHBjC9k2tkXAj4EfAZcDX6ka/jDl2wMehprnSEACEqgfgRIE3FiKslm2dwE/TLIN4cbP9fXDP3hECnhwZp4hAQlIoI4EShDwTcD5qXf7gnRvNwT8VuDkOkIfJSYFPAo9z5WABCRQHwIlCHgF4IEmpG8DDgNWBE4D3gzE8HMRmwIuIo1WQgISkECxzwE/AzgpvaQhesivBb5dQr4VcAlZtA4SkIAEyl6IYzngEOAgYBngM8A7gftzTrwCzjl7xi4BCUjgzwRKGILulc8tgROBjYCrgVel+8W9zqvlfgVcy7QYlAQkIIGBCUyDgAPKSsBHgX2BB4EPAIcPTKsGJyjgGiTBECQgAQmMgcC0CLiBaifgeOBJwCivI/wUEM8fxzPHsdjHFR1y8fqmIfDz0heAxenYbvs6plYBj6HVewkJSEACNSAwbQIO5KsDn05D0cOmYNv0XPGFwPwOAl4fuChNBPsd8HXg7HQvutu+rjEp4GFT5nkSkIAE6kVgGgU8zgzc0EXA7wA2TI9BRZkvBA4Gngd026eAx5khryUBCUigpgQU8GiJ6Sbg6GXfDHw4FbEx8C1gndQD77RPAY+WE8+WgAQkkAUBBTxamiYh4AOA+Fm6zczMzFu0qJh1RUaj79kSkIAEMiaggEdLnkPQo/HzbAlIQAJTS0ABj5b6bgLeAIhJWpsDjUlY5wBHA932OQQ9Wk48WwISkEAWBBTwcGk6FtgZWAu4HbgHeApwHHBW+okrvzE9hhT//b00IavxGFK3fR2jchb0cAnzLAlIQAJ1I6CA65aRHvEo4MwSZrgSkIAEOhBQwJk1DQWcWcIMVwISkIACLqMNKOAy8mgtJCABCdgDzqwNKODMEma4EpCABOwBl9EGFHAZebQWEpCABOwBZ9YGFHBmCTNcCUhAAvaAy2gDCriMPFoLCUhAAvaAM2sDCjizhBmuBCQgAXvAZbSBugv4kk0uWQp6ywVblgHcWkhAAhKYJQL2gGcJ7GxdVgHPFlmvKwEJSGCyBBTwZHmPXJoCHhmhF5CABCRQCwIKuBZp6D+Iugv44idfzAO3PMAmp2/C3Plz+6+YR0pAAhKYMgIKOLOE11nAt555KwtesmAp0TnLz2HjUzdWwpm1L8OVgAQmR0ABT471WEqqq4BDvle//GqWLF7ySD2V8FhS7kUkIIFCCSjgzBJbRwG3k28DqxLOrIEZrgQkMDECCnhiqMdTUN0E3E2+Sng8OfcqEpBAmQQUcGZ5rZuAL1rzIhbfurgnxeXnLs82C7fpeZwHSEACEpgWAgo4s0zXTcD2gDNrQIYrAQnUhoACrk0q+gukbgKOqL0H3F/uPEoCEpBAMwEFnFl7qKOAO0nYCViZNS7DlYAEJkpAAU8U9+iF1VXADQn7HPDoOfYKEpDAdBBQwJnluc4CDpSuhJVZgzJcCUigMgIKuDL0wxVcdwH7NqTh8upZEpDA9BFQwJnlXAFnljDDlYAEiiMwro6GAs6saSjgzBJmuBKQQHEEFPDDKZ1TXGZ7VKjuAp62fFhfCUhg+ggoYAU8fa3eGktAAhKoAQEFrIBr0AwNQQISkMD0EVDACnj6Wr01loAEJFADAgpYAdegGRqCBCQggekjoIAV8PS1emssAQlIoAYEFLACrkEzNAQJSEAC00dAASvg6Wv11lgCEpBADQgoYAVcg2ZoCBKQgASmj4ACVsDT1+qtsQQkIIEaEBjXS2dcinK4ZD4V+CLwBOAuYG9gQculXgvs3/S7tYHzgd2A9YBfAFc27X9p+l3XiFwJa7iEeZYEJCCBcRC49cxbGddrVxXwcBk5DzgROAHYHXgX8Jwel7oKeB9wRhLwFcDqgxavgAcl5vESkIAExkMg5Hv1y69myeIlj1xwzvJz2PjUjZk7f+7AhSjggZGxJnAdsAbwYFqP+hbgeen37a64FfANYB6wWAEPDt0zJCABCVRJoJ18G/EMK2EFPHhGtwC+DDyt6dRLgIOA6Bm32z4H3AO8Pe2MIehrgegFLwucCRwGPNQrHHvAvQi5XwISkMB4CXST7ygSVsCD52lQAa8CRA95a+DqVNyKwAywMPWkTwHOBY5sE84BQPws3WZmZuYtWrRo8Kg9QwISkIAEhiJw0ZoXsfjWGLzsvi0/d3m2WbhNr8Me2a+A+0b1yIGDDkHHBK03Ac/tUtSewF7ALr3CsQfci5D7JSABCYyXgD3g9jyreh/w99IErMYkrBh+fnaHlF+QZkwf17Q/JH5nuh8cveGTgWuAQ3o1GwXci5D7JSABCYyfgPeAH820KgHH/d+Q7+OBu4F45CgeKQrJnpV+Ito47jLgSekecKMG8SjSoeme73Lp3vGBwP29mo0C7kXI/RKQgARmh4CzoP+Sa1UCnp3s9nFVBdwHJA+RgAQkMEsEfA74z2AV8Cw1Mi8rAQlIQALtCbgS1sNcFLCfEAlIQAISmCgB14JWwBNtcBYmAQlIQAIPE1DACtjPggQkIAEJVEBAASvgCpqdRUpAAhKQgAJWwH4KJCABCUigAgIKWAFX0OwsUgISkIAEFLAC9lMgAQlIQAIVEFDACriCZmeREpCABCSggBWwnwIJSEACEqiAgAJWwBU0O4uUgAQkIAEFrID9FEhAAhKQQMYEfB9wZsnzZQyZJcxwJSABCXQgoIAzaxoKOLOEGa4EJCABBVxGG1DAZeTRWkhAAhKwB5xZG1DAmSXMcCUgAQnYAy6jDSjgMvJoLSQgAQnYA86sDSjgzBJmuBKQgATsAZfRBhRwGXm0FhKQgATsAWfWBhRwZgkzXAlIQAL2gMtoAwq4jDxaCwlIQAL2gDNrAwo4s4QZrgQkIAF7wGW0AQVcRh6thQQkIAF7wJm1AQWcWcIMVwISkIA94DLagAIuI4/WQgISkIA94MzagALOLGGGKwEJSMAecBltQAGXkUdrIQEJSMAecGZtQAFnljDDlYAEJGAPuJg2cD9wa41qsyrw+xrF028oOcadY8yRjxzjzjHmXFnnGvc42shcYMV+/2jV7bg5dQtoCuO5EVg7w3rnGHeOMUfTyDHuHGPOlXWucefaRsb251oBjw3l0BfKtRHmGHeOMfvHdeiP1lAn2kaGwjbUSbmyHqqy7U5SwGNDOfSFcm2EOcadY8wKeOiP1lAn2kaGwjbUSbmyHqqyCnhs2MZ6oQOAj4/1ipO5WI5x5xhzZDPHuHOMOVfWucadaxsZ219Ye8BjQ+mFJCABCUhAAv0TUMD9s/JICUhAAhKQwNgIKOCxofRCEpCABCQggf4JKOD+WY165ErAV4CNgfuAhcA+wHXAmsCJwIZAPKe8L3D+qAWO+fzXAl8AXgKcmUHM8Wzgx4DnA38EfgK8Engq8EXgCcBdwN7AgjGzGuVyLwQ+BCwDLAd8NMVbpzbyKWBXYF1gM+CKVOFubOvAvV3c3T6XUa2quXdi3WhjrZ/LOsQcMXSKu9PnMs6pQxsZ5bM78LkKeGBkQ58QH/QdgG8CS4D9gN2B7ZPYfg28H3gO8DVgfWDx0KWN98T1gC8D0V4+kgQcMq5zzP+WBPYvifdawG+B89KXnRMS/3cl5uMlNtzVgu/tqU38FAjuPwNisYFP1oj3tsD1wIXA/CYBd2NbB+7t4u72uYwsVt3OO7GO2Np9LusQc8TQKe5On8s4pw5tZLhP7pBnKeAhwY3htGcDp6cPUayE9ZQkiLj0JcDBwLfHUM6ol4ie2DlAiCp6lJ9IAq5zzKsAt6QFTu5uAhC9mRhxWAN4MH2hiOOel34/KqtRz4/P421plCFGQP4mfWGLL2N31LCN3NAk4G5sIwd14t4cd2vOmj+Xsa8u7bw15k6fyzrFHLE0x93pc9notdepjYz6We7rfAXcF6ZZOeik9Ef1UODmluXUTgW+lb59z0rhA1z0QGA14H3A95KAL6h5zCGus9KQ/05pyD9GFxalnvzTmuofX3YOSt++B8Aya4dGvHGr4g/A44DdgB/XlHfzH9cturCNof4YQakL924Cbnwu9wceXyPurTG3+1zGraE6xdwq4E6fy+8A3dpP9IyL3BRwNWmN3u0uwI7AyjX6kLfS2BT4fBpOiuHwXAS8OXA58Jo03Bz3Kc8Fdk7/XxcRtPKOe74x6nFImgMQtyPii8Sz0vBz85q3dfiSVpqAmz+X99ZMZs2sO30u6y7gTp/LTdJoVZ2+pE3EDAp4Ipj/opD45roHED2d6JHFFr2dmIAV9yhjq8sQdEwSCxnExLDY4j5qDCdGb/iomsYcccYEq98BKwAPpdgvTROajqvxEHQMf8YfoY2aWkzEHWKIP651ayMlDUG3+1zW6bPZzLrb5/Lfa/b3pDnuTp/LdwMx58Eh6Mn7aKpKjJVfXpHke2dTzWNCUDTUxiSs+GMbEyzqMgmrEWqjBxzx1T3muG8d96v/J01oC5E9E/hSir0xCSuGn0N8ddiemP4IbQlck+75xpex6AHHrYq6tZHWYdFoH8G1Hdtu+ybNvjXuTp/LiKsu7bzbsHnz57JOMUcsrXF3+lzelEbYOrWfSbeRiZRnD3gimJcWEm88+k2aPXpPKjZ6llsB8Yc37j3FZJsH0gzp704utL5Lav6g1z3mDYD/SL3hPyWBnZHuQ8aHPO6VRW8+HuO4sm8Cs3/gnqnHGzHHRJsPp15xnXgfm4bzY0QkZm1He45JhDG034ltt32zT/XhEtrFHU8hdPpcxjlVc+/EuplZq4CrjrkT62gjnT6XcU4d2sik2uLSchTwRHFbmAQkIAEJSOBhAgrYliABCUhAAhKogIACrgC6RUpAAhKQgAQUsG1AAhKQgAQkUAEBBVwBdIuUgAQkIAEJKGDbgAQkIAEJSKACAgq4AugWKQEJSEACElDAtgEJSEACEpBABQQUcAXQLbJ4AsumxSmWB36RFgGJJSaH2WJFtF+2nPje9M7gYa43iXPiFYrN621/P71icRJlW4YEsiGggLNJlYFmRCBezXdkWkIylr+MJUXnAbcOUYeGgGO98Hh9ZWyxotc3hrjWpE6J1bv+Kq0d/nxAAU+KvOVkRUABZ5Uug82MQLy9KHqvIaOXpBcqDFqFhoB/ldYHH/T8Ko+PZR5jSVUFXGUWLLu2BBRwbVNjYIUQiHf7/hPwgfSyjUGrpYAHJebxEsiEgALOJFGGmS2BdwFHAF8H5g9Ri34FvA4QveR4sUD8d7stXpywRnoRxR3pgPgbEG/mmkn7/h7YF4h3KMcLIaL3+rbUk4+XQ7wxvcAi3uH6+/S2qXiVX/PbvRpl2wMeIuGeMj0EFPD05NqaTp5AfL7OSxOQhh1C7lfAL033iL8G7NamqvEWmpgQFkPi8d+NLd5Qcy0Qr4OL4eK4Tkg3Xki/A7A68PP01q6zgL9J++P8fwBWAr6azmstVgFPvs1ZYkYEFHBGyTLU7AhEz/HjTVFH77NdT7FbxfoVcEx8incbvwc4vM0FXw6ckiT9sqb9MTwew+SxXQzEvpvT/4eorwJWTpJeAOwN3JX2Pwf4IbAEWBW4r6VcBZxdkzXgSRJQwJOkbVnTROAZwI/S+51DUDHEGz3KQd/z3K+Avw3sCLwAOLsN6I8A7wTenYbEG4c0fn9jGna+reXcS4AQ7RXAc4E/tuwPWccksw3Tu66bdyvgaWrx1nVgAgp4YGSeIIGeBOL53x8AWwCvAV6RhmsPAP6t59l/eUA/Ao7PcdzTjeHiJwBxr7d1+076AhDDxuc27WyI++0tvfXGIdcAT09DzDHU3LxFufcAqwCPTf+tgAdMsIdPLwEFPL25t+azR+BQIBbLiOd2Y0i3MTx8EvDqAYvtR8BPTfdpbwDWb3P9ZkE/Psm6cViI+3HpEae4T928hVhjuDmeY45jWnu/sdhGLLrRel+5cQ17wAMm28Oni4ACnq58W9vZJ7AlcFFadGPTJLuQ8KnAlWkS0yBR9CPgPYFYaSsW6Ni9zcUbE61aBR2yvh74bRpGbj31ecAF6T7v1m2uuxfwpTb3lRXwIBn22KkloICnNvVWfBYIxGSlHwMbAf/YdC827o9el3qSqwH3D1B2PwKOiV4x4avTEpX7A59oI+iQ9WnAfwG7tImpcd7RwFvb7D8KiKHrmPwV95JbN3vAAyTaQ6ePgAKevpxb49kj8GlgP+Az6d/mkmL2c9yjjfvCMTmr360fAZ8P/B8gnsf9WMuFYzWuy4DojR+chsMbhzSGxt+fFgppjSmGzF+ZZj5/sU3AMaEsJBvPDse9ZAXcb1Y9TgKAArYZSGA8BHYCzgH+F9i8zSM58Tzw3wFvAP5jgCJ7CTgWx7g7TYS6NAkxnuGNLR57+mxaiSv+P54PjueEG1vMlo5JWTunBTVaw2pMwIpFN65u2dm6gIcLcQyQVA+VQBBQwLYDCYxOIHq2cX93rfSoTvQ4W7fGcG2n4dxOUfQScMgxntW9BYjh7ZgoFeXHf8eiGXF/NyZQxezoEOqJTY8hxSNHMSnricDClgDi/EVpQY54hCpWxWreGveV4x5yDLG32xyCHr1teYWCCSjggpNr1SZG4OT0qNEhwAc7lNqYKHVhGi7uN7heAo6FMY5Pi2l8Pr2FKYabY/Zy3NuNYefXpvvDsUzlPulZ5HWBmJTVaenKbdOKVzEJK/67dWss4BGTy+K/FXC/GfU4CSQCCtimIIF6E+gl4OhRvyVNhmpedasOtbIHXIcsGENtCSjg2qbGwCSwlECv9wHHUpDx6FP0UqO3WofN9wHXIQvGUHsCCrj2KTLAKSfQEHAzhsbjRrHiVqxEtVxaiaox+apqZLE4RyzS0dh8H3DVGbH8WhJQwLVMi0FJoC8CMdv68iEX+OirAA+SgARmj4ACnj22XlkCEpCABCTQkYACtnFIQAISkIAEKiCggCuAbpESkIAEJCABBWwbkIAEJCABCVRAQAFXAN0iJSABCUhAAgrYNiABCUhAAhKogIACrgC6RUpAAhKQgAQUsG1AAhKQgAQkUAEBBVwBdIuUgAQkIAEJ/H95vXzg/pt3SwAAAABJRU5ErkJggg==\" width=\"599.9999910593034\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# SET UP AXIS\n",
    "waveEx = tauA.waveEx\n",
    "exFlux = tauA.fluxEx\n",
    "exFlux_e = tauA.efluxEx\n",
    "\n",
    "ax3 = plt.figure(figsize=(6,6)).add_subplot(111)\n",
    "ax3.set_title('Excess Flux')\n",
    "ax3.errorbar(waveEx * twinkle._ANG2MICRON, waveEx * exFlux,\n",
    "             yerr=waveEx * exFlux_e, fmt='mD',ms=6)\n",
    "ax3.set_ylabel(r'$\\lambda F_{\\lambda}\\ [erg\\,\\ s^{-1}\\ cm^{-2}] $',fontsize=20, family='sans-serif')\n",
    "ax3.set_xlabel(r'$\\lambda\\ [\\mu m]$', fontsize=20, family='sans-serif')\n",
    "ax3.xaxis.set_major_formatter(mtick.FormatStrFormatter('%d'))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.1 Plot New SED With FAR - IR Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 6. FIT ONE TO TWO BLACKBODIES"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.1 Try fitting 1 blackbody first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# RESET FULL SPECTRUM - bbdust_flux\n",
    "tauA.fullspectrum = tauA.StarPhotosphere[1][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FOR NOW, USING OLD FITTING ROUTINES. WILL EVENTUALLY SWITCH TO PANDAS\n",
    "# OR SCIPY\n",
    "import mosaic_tools as mt\n",
    "ft = mt.FittingTools()\n",
    "\n",
    "# ==============================================\n",
    "# INITIAL GUESS\n",
    "tempdust = 100.\n",
    "SEDx_ang = tauA.StarPhotosphere[0][:]\n",
    "bbdust = {}\n",
    "# ==============================================\n",
    "# CALCULATE BROADBAND PHOTOMETRY ACROSS IR EXCESS WAVELENGTHS\n",
    "for i,band in enumerate(tauA.mags4Dust):\n",
    "    bbdust[band] = stoo.blackbody(SEDx_ang,np.array([tempdust]),1, np.array([band]))[0]\n",
    "\n",
    "bbdust_flux = np.array(zip(*sorted(bbdust.items()))[1]).flatten()\n",
    "\n",
    "FluxNorm_dust = np.average(exFlux / bbdust_flux, weights=1. / exFlux_e)\n",
    "\n",
    "# INITIAL GUESS FOR DUST RADIUS\n",
    "Rad_dust = np.sqrt(FluxNorm_dust / tauA.su2ea_dust)\n",
    "print('Initial guess for dust radius: {} AU'.format(Rad_dust))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "p0_dust = np.array([tempdust, 30])\n",
    "popt, pcov = curve_fit(stoo.blackbody, SEDx_ang, exFlux, p0=p0_dust, sigma=exFlux_e,**)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# FIT BLACKBODY\n",
    "#p0_dust = np.array([tempdust, Rad_dust])\n",
    "p0_dust = np.array([tempdust, Rad_dust])\n",
    "nparams = len(p0_dust)\n",
    "fa_Dust = {'x': SEDx_ang, 'y': exFlux, 'err': exFlux_e,\n",
    "           'func': stoo.blackbody, 'su2ea1': tauA.su2ea_dust, 'bands': tauA.mags4Dust,'bulk':True}\n",
    "\n",
    "#parinfo_dust = [{'value': 0., 'step': 0, 'limits': [0, 0], 'limited': [0, 0], 'fixed': 0} for m in\n",
    "#                range(nparams)]\n",
    "parinfo_dust = [{'value': 0., 'step': 0, 'limits': [0, 0], 'limited': [0, 0], 'fixed': 0} for m in\n",
    "                range(nparams)]\n",
    "\n",
    "for k in range(nparams): parinfo_dust[k]['value'] = p0_dust[k]\n",
    "parinfo_dust[0]['relstep'] = .2\n",
    "parinfo_dust[1]['relstep'] = .2\n",
    "parinfo_dust[0]['limited'] = [1,1]\n",
    "parinfo_dust[0]['limits'] = [10,1000]\n",
    "\n",
    "m_dust = mt.mpfit(ft.deviates_from_model, parinfo=parinfo_dust, functkw=fa_Dust)#,quiet=1)\n",
    "\n",
    "# from scipy.optimize import curve_fit\n",
    "# popt, pcov = \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print out Chi2 and dust temperatures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CALCULATE SINGLE TEMP BLACKBODY FIT\n",
    "p0_dust,errors = m_dust.params, m_dust.perror\n",
    "chi2 = m_dust.fnorm/m_dust.dof\n",
    "print('Chi2={:.2f}, T={:.1f} K , Radius = {:.1f} AU'.format(chi2,*p0_dust)) # , p0_dust, errors \n",
    "#p0_dust[0] = p0_dust[0] * 1.2\n",
    "dust_flux = stoo.blackbody(SEDx_ang, p0_dust,tauA.su2ea_dust)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.2 Plot New SEDs With Fitted Blackbody"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SAVE THE FULL SPECTRUM IN CASE YOU NEED TO RESET\n",
    "fullspectrum_0 = tauA.StarPhotosphere[1][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax5 = plt.figure(figsize=(8,5)).add_subplot(111)\n",
    "\n",
    "tauA.plot_photosphere(ax5, pointsize=8,lcolor='blue',pcolor='blue',lw=0.5)\n",
    "\n",
    "tauA.fullspectrum = tauA.StarPhotosphere[1] + dust_flux\n",
    "\n",
    "plt.plot(SEDx_ang * twinkle._ANG2MICRON, dust_flux * SEDx_ang, 'm-', ls='--',lw=2)\n",
    "\n",
    "ax5.errorbar(waveEx * twinkle._ANG2MICRON, waveEx * exFlux, \n",
    "             yerr=waveEx * exFlux_e, fmt='mD',ms=6)\n",
    "tauA.plot_observedData(ax5,pointsize=8,lw=1.5)\n",
    "ax5.set_title('Excess Flux + BLACKBODY')\n",
    "ax5.set_xlim(.2,300)\n",
    "ax5.set_ylim(1e-15,1e-6)\n",
    "ax5.set_ylabel(r'$\\lambda F_{\\lambda}\\ [erg\\,\\ s^{-1}\\ cm^{-2}] $',fontsize=20, family='sans-serif')\n",
    "ax5.set_xlabel(r'$\\lambda\\ [\\mu m]$', fontsize=20, family='sans-serif')\n",
    "plt.tight_layout()\n",
    "ax5.xaxis.set_major_formatter(mtick.FormatStrFormatter('%d'))\n",
    "ax5.annotate(r'{}'.format(tauA.starname), xy=(0.73, 0.67),xycoords='axes fraction',\n",
    "            fontsize=20, family='Times New Roman')\n",
    "ax5.annotate(r'$T_* = {:.2f}$ K'.format(tauA.StarTemp), xy=(0.73, 0.63),\n",
    "            xycoords='axes fraction',fontsize=10, family='Times New Roman', color='blue')\n",
    "ax5.xaxis.set_major_formatter(mtick.FormatStrFormatter('%d'))\n",
    "plt.loglog()\n",
    "\n",
    "# ===============================================================================================================\n",
    "ax4 = plt.figure(figsize=(8,5)).add_subplot(111)\n",
    "ax4.set_title('Excess Flux + BLACKBODY')\n",
    "ax4.errorbar(waveEx * twinkle._ANG2MICRON, waveEx * exFlux,\n",
    "             yerr=waveEx * exFlux_e, fmt='mD',ms=6)\n",
    "ax4.plot(SEDx_ang * twinkle._ANG2MICRON, dust_flux * SEDx_ang,'m-',lw=2)\n",
    "ax4.set_ylabel(r'$\\lambda F_{\\lambda}\\ [erg\\,\\ s^{-1}\\ cm^{-2}] $',fontsize=20, family='sans-serif')\n",
    "ax4.set_xlabel(r'$\\lambda\\ [\\mu m]$', fontsize=20, family='sans-serif')\n",
    "plt.tight_layout()\n",
    "ax4.xaxis.set_major_formatter(mtick.FormatStrFormatter('%d'))\n",
    "ax4.set_xlim(8,300)\n",
    "ax4.set_ylim(1e-13,1e-9)\n",
    "\n",
    "plt.loglog()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.3 Fitting 2 blackbody's to excess flux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FOR NOW, USING OLD FITTING ROUTINES. WILL EVENTUALLY SWITCH TO PANDAS\n",
    "# OR SCIPY\n",
    "\n",
    "# ==============================================\n",
    "# INITIAL GUESS\n",
    "td1, td2, td3 = 500, 100., 40\n",
    "SEDx_ang = tauA.StarPhotosphere[0]\n",
    "bbdust = {}\n",
    "# ==============================================\n",
    "# CALCULATE BROADBAND PHOTOMETRY ACROSS IR EXCESS WAVELENGTHS\n",
    "for band in tauA.mags4Dust:\n",
    "    bbdust[band] = stoo.NBlackBody(SEDx_ang,np.array([td1,1,td2,1,td3,1]),1, np.array([band]))[0]\n",
    "\n",
    "bbdust_flux = np.array(zip(*sorted(bbdust.items()))[1]).flatten()\n",
    "\n",
    "FluxNorm_dust = np.average(exFlux / bbdust_flux,\n",
    "                           weights=1. / exFlux_e)\n",
    "\n",
    "# INITIAL GUESS FOR DUST RADIUS\n",
    "Rad_dust = np.sqrt(FluxNorm_dust / tauA.su2ea_dust)\n",
    "print(Rad_dust)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIT BLACKBODY\n",
    "p0_dust = np.array([td1,Rad_dust,td2, Rad_dust*2,td3,Rad_dust**3])\n",
    "nparams = len(p0_dust)\n",
    "# fa_Dust = {'x': SEDx_ang, 'y': exFlux, 'err': exFlux_e,\n",
    "#            'func': stoo.doubleBB, 'su2ea1': tauA.su2ea_dust, 'bands': tauA.mags4Dust,'bulk':True}\n",
    "fa_Dust = {'x': SEDx_ang, 'y': exFlux, 'err': exFlux_e,\n",
    "           'func': stoo.NBlackBody, 'su2ea1': tauA.su2ea_dust, 'bands': tauA.mags4Dust,'bulk':True}\n",
    "\n",
    "parinfo_dust = [{'value': 0., 'relstep': 0, 'limits': [0, 0], 'limited': [0, 0], 'fixed': 0} \n",
    "                for m in range(nparams)]\n",
    "\n",
    "for k in range(nparams): parinfo_dust[k]['value'] = p0_dust[k]\n",
    "parinfo_dust[0]['relstep'] = 0.3\n",
    "parinfo_dust[1]['relstep'] = 0.2\n",
    "parinfo_dust[0]['limited'] = [1,1]\n",
    "parinfo_dust[0]['limits'] = [200,1000]\n",
    "\n",
    "parinfo_dust[2]['relstep'] = 0.3\n",
    "parinfo_dust[3]['relstep'] = 0.2\n",
    "parinfo_dust[2]['limited'] = [1,1]\n",
    "parinfo_dust[2]['limits'] = [50,100]\n",
    "\n",
    "parinfo_dust[4]['relstep'] = 0.3\n",
    "parinfo_dust[5]['relstep'] = 0.2\n",
    "parinfo_dust[4]['limited'] = [1,1]\n",
    "parinfo_dust[4]['limits'] = [2,80]\n",
    "\n",
    "\n",
    "m_dust = mt.mpfit(ft.deviates_from_model, parinfo=parinfo_dust, functkw=fa_Dust)#,quiet=1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print out Chi2 and dust temperatures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CALCULATE SINGLE TEMP BLACKBODY FIT\n",
    "p0_dust,errors = m_dust.params, m_dust.perror\n",
    "#chi2 = m_dust.fnorm/m_dust.dof\n",
    "print('Chi2={:.2f}'.format(chi2))# , p0_dust, errors \n",
    "print(p0_dust)\n",
    "dust_flux = stoo.NBlackBody(SEDx_ang, p0_dust,tauA.su2ea_dust)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.2 Plot New SEDs With Fitted Blackbody"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# RESET FULL SPECTRUM - \n",
    "tauA.resetFullSpectrum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax5 = plt.figure(figsize=(8,5)).add_subplot(111)\n",
    "tauA.plot_photosphere(ax5, pointsize=8,lcolor='blue',pcolor='blue',lw=0.5)\n",
    "tauA.fullspectrum = tauA.StarPhotosphere[1] + dust_flux\n",
    "\n",
    "plt.plot(SEDx_ang * twinkle._ANG2MICRON, dust_flux * SEDx_ang, 'm-', ls='--',lw=2)\n",
    "ax5.errorbar(waveEx * twinkle._ANG2MICRON, waveEx * exFlux, \n",
    "             yerr=waveEx * exFlux_e, fmt='mD',ms=6)\n",
    "tauA.plot_observedData(ax5,pointsize=8,lw=1.5)\n",
    "ax5.set_title('Full Spectrum')\n",
    "ax5.set_xlim(.2,300)\n",
    "ax5.set_ylim(1e-15,1e-6)\n",
    "ax5.set_ylabel(r'$\\lambda F_{\\lambda}\\ [erg\\,\\ s^{-1}\\ cm^{-2}] $',fontsize=20, family='sans-serif')\n",
    "ax5.set_xlabel(r'$\\lambda\\ [\\mu m]$', fontsize=20, family='sans-serif')\n",
    "plt.subplots_adjust(left=.19, right=.96, bottom=.15, top=.93)\n",
    "ax5.xaxis.set_major_formatter(mtick.FormatStrFormatter('%d'))\n",
    "ax5.annotate(r'{}'.format(tauA.starname), xy=(0.73, 0.67),xycoords='axes fraction',\n",
    "            fontsize=20, family='Times New Roman')\n",
    "ax5.annotate(r'$T_* = {:.2f}$ K'.format(tauA.StarTemp), xy=(0.73, 0.63),\n",
    "            xycoords='axes fraction',fontsize=10, family='Times New Roman', color='blue')\n",
    "ax5.xaxis.set_major_formatter(mtick.FormatStrFormatter('%d'))\n",
    "plt.loglog()\n",
    "\n",
    "# ===============================================================================================================\n",
    "ax4 = plt.figure(figsize=(8,5)).add_subplot(111)\n",
    "ax4.set_title('Excess Flux')\n",
    "ax4.errorbar(waveEx * twinkle._ANG2MICRON, waveEx * exFlux,\n",
    "             yerr=waveEx * exFlux_e, fmt='mD',ms=6)\n",
    "ax4.plot(SEDx_ang * twinkle._ANG2MICRON, dust_flux * SEDx_ang,'m-',lw=2)\n",
    "ax4.set_ylabel(r'$\\lambda F_{\\lambda}\\ [erg\\,\\ s^{-1}\\ cm^{-2}] $',fontsize=20, family='sans-serif')\n",
    "ax4.set_xlabel(r'$\\lambda\\ [\\mu m]$', fontsize=20, family='sans-serif')\n",
    "plt.subplots_adjust(left=.19, right=.96, bottom=.15, top=.93)\n",
    "ax4.xaxis.set_major_formatter(mtick.FormatStrFormatter('%d'))\n",
    "ax4.set_xlim(8,300)\n",
    "ax4.set_ylim(1e-13,1e-9)\n",
    "\n",
    "plt.loglog()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}