demo_mcfostpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MCFOST Python Library Demo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Assumptions:** This demo assumes you are familiar with the basics of the MCFOST astrophysical radiative transfer code, and have a working copy of that installed on your system. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This demo further assumes you are running this notebook from its location within a copy of mcfost-python git directory; specifically that there's a subdirectory named 'docs/example' which contains the ref2.19.para sample parameter file, and optionally the results of running MCFOST on it to produce an SED and an image at 1.0 microns. You can do that directly from the command line as follows, or we will show below how to invoke this from within Python. \n",
    "\n",
    "   >  mcfost ref2.20.para -rt\n",
    "\n",
    "   >  mcfost ref2.20.para -img 1.0 -rt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%aimport mcfost\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quick test of the extinction calculator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 2.94121201  2.88260425  1.37797835  0.79933598  0.4339158   0.25539264\n",
      "  0.16192279  0.06803995  0.01810171  0.01516129  0.01300451  0.01130458]\n"
     ]
    }
   ],
   "source": [
    "print(mcfost.ccm_extinction(3.1,[0.15,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## General Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This package implements a set of Python object classes to represent the different sorts of datasets one encounters when working with MCFOST: parameter files, model results including SEDs and images, observation results including SEDs and images. The intent is that the model SED class and the observed SED class interfaces look as similar as possible, and likewise for the images. \n",
    "\n",
    "Once you have both model and observed data, it becomes possible to calculate chi^2 for individual fits, and use these to generate the relative conditional probability of a model given the data. This can be used in an MCMC calculation using the emcee package to derive the overall posterior probability distributions. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with Parameter Files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can read in parameter files into a ``Paramfile`` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "par = mcfost.Paramfile('docs/example/ref2.20.para')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting object's attributes provide read/write access to all the settings in that parameter file. There are many, many of these. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distance: 140.0 pc\n",
      "Wavelength range: 0.1  to 3000.0 microns\n",
      "Changing distance to 200 pc\n",
      "Distance: 200 pc\n"
     ]
    }
   ],
   "source": [
    "print(\"Distance:\", par.distance, 'pc')\n",
    "print(\"Wavelength range:\", par.lambda_min, \" to\", par.lambda_max, \"microns\")\n",
    "\n",
    "print(\"Changing distance to 200 pc\")\n",
    "par.distance = 200\n",
    "print(\"Distance:\", par.distance, 'pc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This also includes some quantities which are computed from the parameters, not directly stated in the parameter file.  Currently these computed quantities are only calculated when the Paramfile is read from disk, and don't automatically update if you change some setting. This should be improved eventually. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wavelengths [  1.10859065e-01   1.36242824e-01   1.67438785e-01   2.05777785e-01\n",
      "   2.52895390e-01   3.10801666e-01   3.81966929e-01   4.69427133e-01\n",
      "   5.76913383e-01   7.09011107e-01   8.71355674e-01   1.07087280e+00\n",
      "   1.31607401e+00   1.61741974e+00   1.98776557e+00   2.44291068e+00\n",
      "   3.00227185e+00   3.68971175e+00   4.53455699e+00   5.57284917e+00\n",
      "   6.84888248e+00   8.41709327e+00   1.03443824e+01   1.27129691e+01\n",
      "   1.56238987e+01   1.92013534e+01   2.35979494e+01   2.90012482e+01\n",
      "   3.56417578e+01   4.38027665e+01   5.38324276e+01   6.61586128e+01\n",
      "   8.13071645e+01   9.99243290e+01   1.22804326e+02   1.50923230e+02\n",
      "   1.85480610e+02   2.27950706e+02   2.80145316e+02   3.44291096e+02\n",
      "   4.23124542e+02   5.20008738e+02   6.39076821e+02   7.85408309e+02\n",
      "   9.65245791e+02   1.18626124e+03   1.45788332e+03   1.79169958e+03\n",
      "   2.20195083e+03   2.70613864e+03]\n",
      "Inclinations in RT mode:  [ 45.]\n"
     ]
    }
   ],
   "source": [
    "print(\"Wavelengths\", par.wavelengths)\n",
    "print(\"Inclinations in RT mode: \", par.inclinations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some of the items which can be repeated multiple times in a parameter files are stored as lists of dicts.  For instance, there's a \"density_zones\" list attribute (which can have multiple entries if there are multiple zones defined), each of which in turn contains a \"dust\" list attribute (which can have multiple entries if that zone has more than one type of dust.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[{'dust_mass': 0.001, 'edge': 0.0, 'r_critical': 100.0, 'flaring_exp': 1.125, 'reference_radius': 100.0, 'r_out': 300.0, 'gamma_exp': 0.0, 'zone_type': 1, 'debris_disk_vertical_profile_exponent': 2.0, 'gas_to_dust_ratio': 100.0, 'grain_type': 'Mie', 'dust_nspecies': 1, 'scale_height': 10.0, 'dust': [{'mixing_rule': 2, 'amin': 0.03, 'porosity': 0.0, 'grain_type': 'Mie', 'mass_fraction': 1.0, 'ncomponents': 1, 'filename': 'Draine_Si_sUV.dat', 'volume_fraction': 1.0, 'aexp': 3.5, 'amax': 2000.0, 'ngrains': 100, 'heating': 1}], 'r_in': 1.0, 'surface_density_exp': -0.5}]\n"
     ]
    }
   ],
   "source": [
    "print(par.density_zones)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can all get rather complicated. There are some convenience shorthands defined for accessing the first elements of these lists, for instance ``par.dust`` is equivalent to ``par.density_zones[0]['dust'][0]``. If you want to access more than the first one you must use the full syntax. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cautionary Note:** The API for this all could use some cleanup; the names used for the various different parameter settings are not entirely consistent. We should think some about how much we can try for consistent naming conventions between the Python, IDL, Yorick variable names and the name/descriptions used in the software manual.  And furthermore the Python code here is not yet consistent in terms of using objects with attributes versus dicts. This needs to be made better, time permitting. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can view a parameter file in text form just by casting it to a string, and there is a \".writeto()\" method to save it to a file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.20                      mcfost version\n",
      "\n",
      "#Number of photon packages\n",
      "  1.28e+05                nbr_photons_eq_th  : T computation\n",
      "  1280                    nbr_photons_lambda : SED computation\n",
      "  1.28e+05                nbr_photons_image : images computation\n",
      "\n",
      "#Wavelength\n",
      "  50  0.1   3000          n_lambda, lambda_min, lambda_max [microns]\n",
      "  T T T                   compute temperature?, compute sed?, use default wavelength grid ?\n",
      "  IMLup.lambda           wavelength file (if previous parameter is F)\n",
      "  F T                     separation of different contributions?, stokes parameters?\n",
      "\n",
      "#Grid geometry and size\n",
      "  1                       1 = cylindrical, 2 = spherical, 3 = Voronoi tesselation (this is in beta, please ask Christophe)\n",
      "  100  70  1 20           n_rad (log distribution), nz (or n_theta), n_az, n_rad_in\n",
      "\n",
      "#Maps\n",
      "  101 101 3000.0        grid (nx,ny), size [AU]\n",
      "  10 1               MC : N_bin_incl, N_bin_az\n",
      "  45.0  45.0   1 F    RT: imin, imax, n_incl, centered ?\n",
      "  0.0   0.0    1    RT: az_min, az_max, n_az\n",
      "  200.00                 distance (pc)\n",
      "  0.00                    disk PA\n",
      "  \n",
      "#Scattering method\n",
      "  0                      0=auto, 1=grain prop, 2=cell prop\n",
      "  1                      1=Mie, 2=hg (2 implies the loss of polarizarion)\n",
      "\n",
      "#Symmetries\n",
      "  T                      image symmetry\n",
      "  T                      central symmetry\n",
      "  T                      axial symmetry (important only if N_phi > 1)\n",
      "\n",
      "#Disk physics\n",
      "  0  0.50    1.00        dust_settling (0=no settling, 1=parametric, 2=Dubrulle, 3=Fromang), exp_strat, a_strat (for parametric settling)\n",
      "  F                       dust radial migration\n",
      "  F                      sublimate dust\n",
      "  F                       hydrostatic equilibrium\n",
      "  F  1e-05                 viscous heating, alpha_viscosity\n",
      "\n",
      "#Number of zones : 1 zone = 1 density structure + corresponding grain properties\n",
      "  1\n",
      "      \n",
      "#Density structure\n",
      "  1                       zone type : 1 = disk, 2 = tapered-edge disk, 3 = envelope, 4 = debris disk, 5 = wall\n",
      "  1.00e-03   100.0        dust mass,  gas-to-dust mass ratio\n",
      "  10.0   100.0  2.0             scale height, reference radius (AU), unused for envelope, vertical profile exponent (only for debris disk)\n",
      "  1.0     0.0    300.0  100.0   Rin, edge, Rout, Rc (AU) Rc is only used for tappered-edge & debris disks (Rout set to 8*Rc if Rout==0)\n",
      "  1.125                     flaring exponent, unused for envelope\n",
      "  -0.500 0.000                 surface density exponent (or -gamma for tappered-edge disk or volume density for envelope), usually < 0, -gamma_exp (or alpha_in & alpha_out for debris disk)\n",
      "        \n",
      "#Cavity : everything is empty above the surface\n",
      "  F                        cavity ?\n",
      "  15.0   50.0              height, reference radius (AU)\n",
      "  1.5                      flaring exponent\n",
      "        \n",
      "#Grain properties\n",
      "  1                       Number of species\n",
      "  Mie 1  2 0.00  1.00  0.9     Grain type (Mie or DHS), N_components, mixing rule (1 = EMT or 2 = coating),  porosity, mass fraction, Vmax (for DHS)\n",
      "  Draine_Si_sUV.dat  1.0    Optical indices file, volume fraction\n",
      "  1                       Heating method : 1 = RE + LTE, 2 = RE + NLTE, 3 = NRE\n",
      "  0.030  2000.0   3.5 100   amin, amax, aexp, nbr_grains\n",
      "      \n",
      "#Molecular RT settings\n",
      "  T T T 15.              lpop, laccurate_pop, LTE, profile width\n",
      "  0.2                    v_turb (delta)\n",
      "  1                      nmol\n",
      "  co@xpol.dat 6          molecular data filename, level_max\n",
      "  1.0 50                 vmax (km.s-1), n_speed\n",
      "  T 1.e-6 abundance.fits.gz   cst molecule abundance ?, abundance, abundance file\n",
      "  T  3                   ray tracing ?,  number of lines in ray-tracing\n",
      "  1 2 3                  transition numbers\n",
      "      \n",
      "#Star properties\n",
      "  1 Number of stars\n",
      "  4000.0    2.0   1.0   0.0   0.0   0.0 T       Temp, radius (solar radius),M (solar mass),x,y,z (AU), is a blackbody?\n",
      "  b'lte4000-3.5.NextGen.fits.gz'\n",
      "  0.1    2.2       fUV, slope_fUV \n"
     ]
    }
   ],
   "source": [
    "print(par)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running MCFOST calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MCFOST calculations may be executed for a given parameter file or for many parameter files at once. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that when you run a parameter file, the wrapper function may optionally move that parameter file to a new subdirectory before running it. This is useful if you're running a grid or MCMC and want each parameter file's results in its own subdirectory (the typical case). But you can disable this behavior if it's not desired. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, the following line will create an SED and two images for the example reference parameter file included with MCFOST and with this Python package. Warning, this will take a long time to run:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:mcfost:Running MCFOST for: docs/example/ref2.20.para\n",
      "INFO:mcfost:Computing SED for docs/example/ref2.20.para\n",
      "INFO:mcfost:SED results written to /Users/cpinte/mcfost-python/docs/example/data_th\n",
      "INFO:mcfost:Computing image at 0.8 microns for docs/example/ref2.20.para\n",
      "INFO:mcfost:Computing image at 1.0 microns for docs/example/ref2.20.para\n",
      "INFO:mcfost:Calculation complete.\n"
     ]
    }
   ],
   "source": [
    "type('docs/example/ref2.20.para')\n",
    "res = mcfost.run_one_file('docs/example/ref2.20.para', wavelengths=[0.8, 1.0], move_to_subdir=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with Model Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ModelResults class provides an interface to work with MCFOST model results for *a single model*. You instantiate it with a directory path to an MCFOST results directory (i.e. the parent directory that contains 'data_th' as a subdirectory.)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<MCFOST ModelResults in directory '/Users/cpinte/mcfost-python/docs/example'>\n",
      "['0.8', '1.0']\n",
      "Model results in /Users/cpinte/mcfost-python/docs/example for ref2.20.para\n",
      "    Model has 1 inclinations from 45.0 to 45.0\n",
      "    SED computed from 0.1 - 3000.0 microns using 50 wavelengths\n",
      "    Images computed for 2 wavelengths: [ 0.8  1. ] micron\n"
     ]
    }
   ],
   "source": [
    "model = mcfost.ModelResults('docs/example/')\n",
    "print(model)\n",
    "print(list(model._wavelengths_lookup.values()))\n",
    "model.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with Model Results: SEDs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model contains an ``sed`` attribute which is an object representing the calculated SED, and an ``images`` attribute which is a dict of objects representing all available image wavelengths. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the ``sed`` object, there are attributes ``.wavelength``, ``.flux``, ``.nu_fnu``, ``.uncert`` and so forth. The results are returned as ``astropy.Quantity`` arrays with the appropriate units.  ``wavelength`` and ``inclination`` are 1D arrays, while ``flux`` and so forth are 2D arrays with shape equal to (# of inclinations, # of wavelengths)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x10bb0a7b8>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAEiCAYAAAALRPCFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecFdX9//HXh6qCCthBBEEs2I0CVlZsoEaJioho7Ikx\nlhg1dlnjV8Vo1DRLbCBWLLE3fsoasaHRKMZeUIpiAaSpIPv5/XHmssN17+7d3bk79959Px+Pfeyd\nuXNnPjN37nzmnDlzxtwdERERKQ+t0g5AREREkqPELiIiUkaU2EVERMqIEruIiEgZUWIXEREpI0rs\nIiIiZSSVxG5mPc2s2szqXb6ZHWlmzzVHXPkwsx3N7AMzm29m+6UdTzkxs53N7N2042goM9vIzP5r\nZvPM7MS048lXmr8tMxtjZhelsexiFR0TexXbvJoQQ97H+YSXW2lm4wo07xXN7GEzm2tm483sUDN7\nshDLaop8EutUM/vBzFbLGv969KWtV7jw6mdmx5jZO9FB9Qsze9TMOkbvjYlinx/7ez16L7PTZcZ/\nEX1hu9ezyD8Cf3X3ld39oUKvX23MbISZ3W5mA81sWi3vV5nZMSnEdb2ZHdfYz7v7c+6+cZ7Lqqht\n3fP43Agzu73h0dXpD8DT7r6Ku/+9sTMp5ME4jYOsmXWt4zvy6K+Qyy/Ed100zOxsM7u4ifNIbL+I\ncsWgRnzuyTyOuw3RqP0qz21xELAm0MXdD3b3O9x9r9g8Uj+hgvxK7A58DIzIjDCzzYEVKfAPsz5m\nNhC4GDjE3VcBNgHuik3iwGVREs78bZ01m1XdfWVgC2AC8C8zO6KOxa4HvJ0jHjMza+z6NMA+wKNA\nrmUletBswHoNjuIqZpltl6Qe5NgnGqHQ+09z7J8ZewOP1/F+oWMpxHddTPYmufVL4rvwhs7HzDoA\nPwOeTWD5jWZmbeKDdUzaA3jf3avrml0yUTWBu9f5B3wCnAtMjo27AjgHqAbWi8atCtwKfAlMjT5j\n0Xutos98BXwE/Db6bKvYZ28CZgLTgYti7x0JPJcjttOBf9UR+y3AH3O81zMeQ2z8acAXOT7zEbAU\nWATMA9oBVcD/Ac9H43sBOwCvAHOBycD2sXlURev3PDAfeAhYDbgd+Daavkcd69QK+ALoAlQA02qZ\nZiJwdPS6H/BqNO8vgD/HphsAvADMAf4LDMyKM75evaPv4qNo3T8GDo1NvwXwRmz4OEKymwf8D9gq\nGj8VOCsaNxu4GWgfvbfc+kTTnga8EW3Lu4D2QAfgu+i7mB8tY23CD+os4EPga+BuoHOObZf5/o8E\nPoti+TWwHfBmtE3+lrVds9dpa+AZ4MconnnABtT9W9iAcBCbS/g93BmN/3cUz4JonYbl+P7r2n5v\nAfvGpm0bbYetonWsjm2vAcARwHPA5dG8PgYGxz7flbB/fgN8ABwbe68SGA+Mjeb3FvCzrFjvB4ZG\nr7cGXoumvQu4E7goa9t+EC3rQWCd2HubEk66v4m+v7Pz2LeXfdd17euE3+pXwLrR8JbRttgwGs7s\nT5nvfGhsGUcSfh9XRvP9CNiemn1qFvDL2PRjgOuAp6L5VREdP6P3q4Fe0ev2hGPmp9F6XAusEJu2\nczT/zH51BjXHz6Oz5rUP8Hq0nT4DRsXmE98v5gP9Cb/1Zwj7zlfAbYQCUF15Yhw1x8b5hGNzz2je\nv4zW4yvgnKzP7Qc80MRj1fqE39S8aNv+DRhXT7yZ2I6OYns2+r/ctsj6zIXAD8Di6P2jieUnavkN\nA6sDj0RxfxNNYzliymynr4HzCL/1QdF7c2NxLSCWe2udV10rH83wE2A34F1gY6A1MI1Qco0n9luB\nfxEOuj2A96hJLscD7wDdCDvkxGgnyCTvfxF23BWBNYCXgV/Ffjy5EvtO0Y5UCexIdICLvX8LsYNH\nji82O7H3isZvVMf2GBQbroq+gE0IB5O1oi9xZDR8COFA0Tk2/fuEnXEVwsHiPWBQtG3HAjfX8X0M\nAF6IXldQf2J/ERgZvV6JaGeNvouviQ7kwO7R8Go51mtVwg+uT/T+WkDf2DLPAi6OXg8jHGB+Fg33\nju0nUwmJM7MvTMp8R9nrE23rlwhJuzMhqf46em9g9roDpxB+/F0JSe064I4c2y7z/V9DOEHbA/ie\nkIxWj+YxC9glj3Vatr3z+C3cSU1iagfsUNuBvY7vv67tdwZwV2za/YlOtqI4ltvfCb+txcAxhJOi\n44EZsff/Dfw9inNLwonKrtF7lYSTmcHRZy8BXox9ti3hQN4h+vyn0ffTGjgwWu4fo2kHRdNuFU37\nV+DZ6L2Vgc+BU6P3OgLb1bVv1/Jd17ev/x/wNOH4MwU4ITafg4C1o9cHEw6qa8W23xLCCZIRTtg/\nIySWtoR9ah6wUjT9mGh4p2hdriZ2bGP5ZHwV8ADQKVrnh4BLYtMeAtwevR5MSIR9o+1wR9a8BgKb\nRq83j6bdv479ojfhmN+W8Ft4Frgqz1wRPzb2jOZ9PeFEZQvCb2zj2DTXAcc18Vj1IuEkqC2wc7SN\nb60n1kxsY6LvvX1t26KWz42Kz5us/ETWbxi4lJDbWkd/O+aYb19C0t4hWo/LCb+RQbVMewnh+Nw6\nZ5x5flm7EUodl0Q70ZNRkNWEBN+acCYT/8J+BUyMXj9DlKij4T0yG5CQIL5n+bPREcAztW24WuIb\nTNjp50Qb5s/UnDCMIRx85sT+bsn6YrMT+wrR+O3r2B7xnXciUBkbPhx4KeszLwBHxKY/O/beFcCj\nseF9gdfrWN+LgHOj1xXUn9ifJRyEV8+a5kyydn7gCaISRi3r1SHafgcAK9ayzH9ndtpo/zipju0X\n3xeGAB/Wtj7RtPFagcuAa3OtOyHxx7+bdQg/jsz+EN92me8/XjL8mlhJGbgXODmPdZoIHBO9ru+3\nMJZwoOtWy3zySex1bb+uhINax1j8p+fa3wm/rQ9iwytF06wJdCfURHSIvX8JNb+fSuCp2Ht9gUWx\n4d2ACdHrXYidMETjnqcmsd8EjM7a1xYTDrQjgP/k2Ba17tu1fNf17ettCCXFKcBj9Wz/14H9Ytvv\n/dh7m0fbb42sfWqL6PUYlj/R7BBt427x759wkrCA5RPE9sDHseFx1CTBm1k+6fepa18inFBcmWu/\nqGX6ocBrdW2X2L5ZW2LvGhv3MjA8NvxpbP0bfKwi5J8lxI5JhNrPfEvsPWsZV9e2qIzPm/oT+4WE\nE7Te9cRzAdGJWjS8IuE4MihruuHRdl6trvnl22DCMzsS4ez0Vpa/jrA64Szj09i4zwhnWhAOsNOy\n3svoEX32czObY2ZzCGdxa+QVmPsT7r6fu3cmlFCOBI6NxX25u3eO/R1VzywzMc/OZ/mR+Lp1Zfn1\ng7BdusaGZ8Vef08oCcWHO9axrCHAY9HrHwnbLltbws4OoTS2IfCOmU02s32i8T2AYZltHm33HQml\n45+sl7svJOxUxwMzzewRM9sIwMw6EWpzXogmX5dQLZlL9r7QNdeEhNJFxnfUvW16EtpIZNbnbcI2\nWit6P77tMuLfxXe1DGeWV986efS/vt/CHwi/nclm9paZ5dwfzezxWOPOEbG3at1+7j6TkDAPir6T\nwYSDXF2WbV93XxS97BjNc3b0vde2HrD8tloErBBreLQ3Ndu6KzAja7nx7bNOfDha5jfRstYlXCKo\nTa59G5b/ruvc1939R8IJ16aEgsEyZvbLqKFw5nObES6d1bYNvovm91XWuMw+5IRan/h6zuan+/8a\nhJOs/8SW+zhh3yLaxrsTkhvUfXzFzPqb2UQz+9LM5hIuOS3XGDpr+rXM7C4zm25m3xKO/Tmnz0P8\nN7yIcEKTaav1rbtn9o3GHKu6AnPc/bvYMuL7Vn0a3AC3gS4nXMp5ysw+MrMzc0y3DsvvG98RfgPL\nmNnWhNqgoe7+DXVoU9ebce7+mZl9TPjBHJ319teERNKTUOUO4UwqE+jn0TCx9zKmEc5MVvO6GyTk\nE+MzZvYM4QfaWL8AZrn7ew1ZdOz1DEKpNq4HuRsReY7xP2FmaxNKmK9Hoz4DVjezDpkDcNTIrQfR\nzu3uHwKHRu8dCNxr4Q6Hzwhnnr/Kc71w96cIO2h7QqPFGwilsb0IrcIz008jXEvOJXtfmFnniucR\nW+Qz4Ch3fzH7jVq2XUPVt04Zdf4W3H0WoQSPme0I/D8ze9bdf5K83H1IjmXUtf3GEk5s2xKqoj/P\nzC6P2ONmAl3MrKO7L8hejzwMIfyWIPz+u2W934NwwMssq2fmjahB1WrRsqYRqp1/Ise+3YVw2Sj7\nd5JzXzezboQS083AlWa2nbsvNrMewD8JlwpedHeP7qppbOMoI9SEZJbbkdDeI3v//5pwQtA39v3F\nbQd8Gju413V8hVA1/1dgr2i9riI6SaD2/eISwqXSzdx9rpkNJSSU+jR0H1uu8V9jjlXRd9TZzFaK\nnZj2iOLPh+d4nc/09U8cfjunA6eb2abAM2b2irs/kzXp58BGmQEzW5HYyZSZrUm4vHeCu79R33Ib\neovDMYSqgfjZEe6+lNCQ5mIz6xht7FMJjS6I3jvZzLqZWWfC9djMZz8nNHi40sxWNrNWZtbbzHap\nLxgz28/MhptZ56jhdj/C9aSXMpNQ/4/QonmtZeEe5AuAs+tbdm3ziDwGbGjhVps2ZjacUJp9JMf0\nDTlIDCF2guDunxGqti4zsw5Rwj2DUI35EoCZHWZmmdqPbwk75lLCd/NzM9vTzFqb2QoWbiGLH4CX\nxWZma5rZ/tFBdwmwkJofT3br3BsJO/I20feygdXcFmnACdG+0IVwiSd+J0O+ZgGrmdkqsXHXAZdk\nlmVma1hNXwPLbbsGyGyDutZp2XT1/RbMbJiZrRt9Zi7h+8ic0M4iXN+sL566tt8DwDbAyYSatYyv\nouXUN3+i9ZhGqIG51Mzam9kWhBP62+r+JJjZ+oT2LpmT4xeAH83sZDNra2YHEJJTxp3AUWa2ZbQP\nX0K4nPUZYb9ax8xOieJYOfqd59q3q/npd51zX49OhMcAN7r7sYQDbOb++g7RPL8GWlmoXdksn+1X\nh70t9IXRLlrOi7ESKwBRAecG4OrM+kWx7pmZB8sfT8YDR5rZJma2EuE6cFxHQql2cbTtDqUmQdW2\nX3Qk/L7nRceDM/Jct3z237ghxI4bjTlWufunhMsoF0b71k6Ey5mNkc9vpL7j9XLbwMz2iY4VRrhM\ntpTaTzruJazj9tG+UUlNbmoTvX+bu9+bz4o0KLG7+8fu/lp8VOz1SYSd4WNCS9vbCY3XIOykTxJa\nN78K3Jf12V8SGpO8TaiauoeaKmEn91nSHEJr2vcJO8I44E/ufmfss3+w5e9j/zJrHnPNbAGhQdJg\n4CB3H1PXdqjFsvjcfTZhxzqNcEA4ndBSeXZt0+dYv1zruw8/rUoeTrgm+iGhhLMrsI+7L47e3wt4\ny8zmExrkHOLuP7j7dMKli3MIlwI+i2KO77jxOFoREtQMQhXRzsBvoh12T2qqBYl2vosJJYV5hAZp\nnWPzvINwMvcRoSX0/+Wx7pn3PFrGu4SE8LGZzbZQIv8Lob3FU2Y2j9Copl8d2y7vM/R61il7Xj/5\nLbj7zdF72wIvRd/Hg4Rr+FOj9yqBsRaqGw+qI56c28/dvyf8vnpGMWbGL4rifz7aXv2pf98bEc1n\nZjSvC2Iljbo+u9xtZu6+hFCLdSRh3zk4ijHz/tPA+dG4mYSGpYdE780ntMn5OSHpvk9oXwE59m2y\nvus69vVWhBOg1aPlAxxFOMnY0d3fJlTNv0ioTt6M0Fgxvr75/nYz791BSLzfEO4UOCzHZ88k/KZf\nslAdPoFQTQ3LX+bA3Z8gXDd/Jto+T2fN6wTgj9Fv4nzC3SKZz2bvF/0I14W3IRxTH+anx+tcLgXO\ni/bf39eyTstYuFTUl5rLd9DwY1Umfx1KaM0/m1AwG5tHrD+JrbZtYaHTrPlZn6nr+F1JzW94GKG9\nwwRC+68XgH+4+7PRNnjMzM6Klv024bhxF+E3MD9a1x8Il6N2An4Xy2PzYgWEn8jcKiFZolLHuYTb\nPIZlD6cQTxvCgW39WNVo6qIDwV/dfUCe039CaGiWXRVVMMW67Rojn+1nZhcAG7j7L5svsuWW/yjh\nVsEn6p04+WUX7XdtZrcA0939/Honzj2PtQgN2bIvbZQUMzsYOMDda73M0tJZuEwzh/A7bkibAUB9\nxefk7p9EVXO1DqegM3BesR2sCGer2VV/xaZYt13iour5ownXhtNSFf2loZi/6yQ6LlkF+H29UxW/\nOYRSuUTM7OdmtpKFy51XAG82JqlDC0jsZnazmc0ysylZ4web2bsW+n3P1VKxaLj7V+5+fdpxZHP3\nV9y96PpKjivWbZc0C935fka4ZWtSfdMXirtfHl0SSGPZxfxd13VZMb8ZuH/g7nfXP2XyzGy9rMua\neVUL18bdJ7j7y4WKFcDMRuaId0r9n07FfoRLnTMI1+kbXZtR9lXxZrYz4Z7QW91982hca0KnIbsT\nNuIrhOuJ2xKuLV3u4dYhzOyeeNV79rCIiEgxKfsSu7s/R6j2ietH6NRjatSw5y5CT0zj3P1Ud59p\nZl3M7DpgKzM7M3u4mVdDREQkL3nfx15murF8xwTTCa0ql4lasR+f9bns4Z8ws/KuAhERKRB3T/8B\nKmWg7EvsORQ0+Y4aNYqJEyfW2f1ivn+jRo1KdPq63q/tvexxDRluaOzaFtoW2hYtb1tMnDiRUaOK\nvf1taWmpJfYZxHqAil7n26NWvSorK5OaFRUVFYlOX9f7tb2XPa6hw0nStmj8vLUt8p9e2yK/95Pa\nFpm/Cy+8sM5YpAGSPFss1j9CJxtTYsNtCJ179CR0jPNfYJOEluWjRo3yiRMneks3atSotEMoGtoW\nNbQtamhbuE+cONFHjRrlIR2lny/K4a8ltIq/k9DN7GqEnnwucPdbzGwIobem1sBN7n5pQsvzct+m\n+aqqqipo6aSUaFvU0LaooW1Rw8xwXWNPRNkn9uZmZj5q1Khl1UsiIpJbVVUVVVVVXHjhhUrsCVFi\nT5hK7CKlKzz6QAqttmOkSuzJaamN50REaqUT88LSyVPhtdTb3QqqsrKSqqqqtMMQESl6VVVVid5J\nJKqKT5yq4kVKV1QdnHYYZS3XNlZVfHJUYhcRESkjSuwFoKp4EZH8qCo+eaqKT5iq4kVKV0uuin/2\n2Wfp168fZsYrr7zCzjvv3OB5PPjggyxZsoQ999yTVVZZpdZpVBVfeCqxi4i0QK+99hqjR49eNnzE\nEUfQoUMHevTowezZs2v9zAMPPMAll1zC6NGjGTdu3HLvPfTQQwwaNIiDDjqIBx54oKCxS910u1sB\nVFZWqoMaESmYs846i0GDBrHnnns26vPV1dWcf/759O9f81DLc889l8GDB9O1a1dat279k898++23\nXHTRRfznP/8BYPvtt2fIkCGsvvrqy6bJlMQbcktbpoMaSY6q4hOmqniR0tVSquLvuecePv30UxYu\nXLjsyWpjx47liCOOyPmZhx9+mPHjxy8rqR9//PHstttuDBs2bNk0DzzwAEuWLGGPPfagU6dOtc5H\nVfGFpxK7iEiJ+/jjj7nhhhtyvj9gwAD2339/AL7++mtatWrFGmuswcKFC5dNM3nyZKqrq5k9ezZ9\n+vRhv/32W24e06dPXy5Zd+rUiQ8++GC5aYYOHZrE6kgTKbGLiDTArrs2fR4TJzbuc7Nnz2bChAnc\nd999jB8/ftn4Xr16ceml+T3H6v777+e4447j1ltvXW78McccwzbbbAPAVlttxS677LJcIp87dy4r\nrLDCsuF27dqxYMGCxq2IFJQSu4hIAzQ2KSfhtddeY6+99uKKK65o1OcnT55M//79a60O32qrrZa9\n7ty5M1VVVcuVwFdeeWW++eabZcPfffcda621VqPikMJSYi8ANZ4TkULYfffdufrqqznyyCOXG59v\nVfzLL7/MokWLeOKJJ3j++ef57rvvePDBB5k/fz6PPfYYd9xxBwALFiygTZvl00Pv3r159dVXlw1/\n/fXXy0r4TaHGc8lT47mEqfGcSOkqhcZz/fv356mnnmLSpEnss88+jZ5PZWUlZsaoUaOYNGkSS5cu\nZeDAgSxcuJBNN92Ut99+m5VWWomPPvqIXr16sWjRIgYMGMCUKVMA2HLLLZkwYQJrrrlmg5arxnOF\np8SeMCV2kdJVCol9xIgR7Lvvvuy5556sscYajZrH+PHjGT16NGbG2WefzUEHHcTtt9/OV199xSef\nfMKhhx667Fa4bbbZhptuuomtt96acePG8emnn1JdXU3v3r0ZOXJkg5etxF54SuwJU2IXKV2lkNhL\nnRJ74annORERkTKixC4iIlJGlNhFRETKiBJ7AeixrSIi+dFjW5OnxnMJU+M5kdKlxnOFp8ZzhacS\nu4iISBlRz3MiIjENeeSoSDFSYhcRiagaXsqBquJFRETKiBJ7Dma2vpndaGb3RMP7m9k/zewuM9sj\n7fhERERqo1bx9TCze9x9WGy4E3CFux+bY/oW1yr+yy/hxRfhu++gQwdYaSXo2DG87tgRuneH1q3T\njlJEiplaxSen7K+xm9nNwD7Al+6+eWz8YOBqoDVwo7tflucszwP+nnigJcQdpk6FSZPg+efhiy9g\nwADo1AlmzICFC2HBgvD/229h8WLYf3/YZ58wjYiIFE7Zl9jNbGdgAXBrJrGbWWvgPWB3YAbwCjAC\n2BbYBrjc3WdG097j7sMsNJUdDTzl7k/XsbyyLbFXV8P48fDII/Djj7DTTrDjjrDFFnWXyD/4AP71\nL3juufCZX/wCNtyw+eIWkeKnEntyyj6xA5hZT+DhWGLfHhjl7oOj4bMA3H107DNdgEuA3YCbgIXA\nEYSTgP+6+/U5llWWiX3RIhg9GubMgZNOgj59oKF3BX37LTz6KDz4IKy5Jhx3XDgpEBFRYk9O2VfF\n59ANmBYbng70j0/g7rOB47M+97d8Zh7vHrGiooKKiorGxFg0Zs6E886Dvn3h/POhbdvGzWfVVeHQ\nQ2H4cKiqggsvhD32gKOPhnbtEg1ZRIpcVVWVut4ukJaa2AtepC6HhA7w+uvwf/8Hhx0GQ4c2vJRe\nm9atYbfd4Gc/gyuvhOOPh7PPDrUAItIyZI6RSvDJa6lV8QOAylhV/NlAdQMa0NW1rLKoincP18Vv\nvz2U1rfeunDLmTABrr0WDjgglOjVgl6k5VFVfHJa6n3srwJ9zKynmbUDhgMPJTXzUn+626JFcNll\noZHc3/5WuKQOoQZgzz3h+uvhjTfC9ftZswq3PBEpLnq6W/LKvsRuZncCA4HVgC+BC9z9FjMbQs3t\nbje5+6UJLa+kS+xTpsCll4Zk/tvfhnvSm0um1f2DD8Kf/wxduzbfskUkXSqxJ6fsE3tzK9XEvmQJ\njBkDTzwBv/99uI0tLQ88AHfdFa6/K7mLtAxK7MlpqY3nCqqysrKkGs9NnQoXXxxuQbvxRujcOd14\nhg6FVq3g1FNDyX3dddONR0QKR43nkqcSe8JKqcReXQ333RcayB13HOy9dzKt3pPyyCNw660huXfv\nnnY0IlJIKrEnRyX2AiiFEvv06fCnP4VEfs01xVnlve++oYX8aafBFVfAeuulHZGIJE0l9uSpxJ6w\nYi+xV1eH29jGjYPDDw/du7Yq8nsjnnwyXCL485+V3EXKlUrsyVFiT1gxJ/YZM0Ip3R3+8IfSunb9\n+OMwdiz85S+w1lppRyMiSVNiT46q4gug2Kriq6tDS/Nbbw09yB1wQPGX0rMNGRKeFnfGGSG5p93A\nT0SSoar45KnEnrBiK7HPmAGXXx6S+xlnlH4jtFtugRdegKuvDs97F5HyoBJ7cpTYE1YsiT1+Lb1U\nS+m1cQ+94X30Uegdb4UV0o5IRJKgxJ4cJfaEFUNiz7R4h9K7lp6P6urQO96CBXDRRdBGF5RESp4S\ne3LKoAxXfNLqK766Gu69F048EQYODNXV5ZbUIdQ8nHlmuFVv9Oiw3iJSmtRXfPJUYk9YWiX2mTND\n1XQptnhvrB9+CAl+o43gN79JOxoRaQqV2JOjxJ6w5k7s1dXw8MOhUdnIkXDggeVxLT1f8+bBCSeE\nx73uvXfa0YhIYymxJ0dXJ0vYF1+EFu/ffx8alJV6i/fGWGWVcL39lFNC73lbbZV2RCIi6WpBZbvm\nU+hr7O6hH/Xjj4fttmu5ST2je3c47zz44x9Dw0ERKR26xp48VcUnrNBV8d98E/pNnz0bzj4bevYs\n2KJKzsMPwz33wD/+ASuvnHY0Nd58M9ylMG8erL46rLZazd/qq8P228Paa6cdpUi6VBWfHCX2hBUy\nsVdVwV//Cj//ebg3vW3bgiympP3jH/Dxx6EhYdq3wS1ZEto+PPVUeMb9ZpvB11+Hk7Ovvw4nZ59/\nDs89B/37w4gRsP766cYskhYl9uQosSesEIl9/vyQ0N99F845BzbZJNHZl5Xqajj33PBs+d/9Lr3H\n0Gaecb/WWuHpdHV1gbtgATz0UHiE7kYbhYaAm23WbKGKFAUl9uQosScs6cT+2mvhXu2ddoJf/Uo9\nreVj4UI46STYc0845JDmXXZ1Ndx/P9x2W/i+hgzJ/+Tihx/Ck+zuvjtU02+4IbRrF2pm2rULf+3b\nwy67QKdOhV0PkeamxJ4cJfaEJZXYq6vhjjvCw1vOOgu23TaB4FqQr76Ck0+GX/4yJNfmMGMGXHll\nSNDnnNP4Z9wvXRr6w//yS1i8uOZvyZJQjT9lSmgsuMUWycYvkiYl9uQosScsicS+YEG4hevbb6Gy\nMjSwkoabNi1Ux//+97DjjoVbzpIlMH58aLg3cmTol79168It7+WXQ2O8Aw4I1+VbUr8FUr6U2JOj\nQ0IBNOV2t48+Crexrb02XHWVknpTdO8Ol1wS7iJ4443CLOOtt+DXvw7/r7sOhg0rbFKH0NDuuutC\ngj/rLJg7t7DLEykk3e6WPJXYE9aUEvuECaFV929/C3vskXBgLdhrr4WHxfzpT9CnTzLzXLgQbrgB\nJk2q6Zu/uRvqLV1a0+r+3HNhyy2bd/kiSVKJPTlK7AlrTGL/7rtQAvvPf+DCC6F37wIF14L9+9/h\nzoK//AVP5rqyAAAfbElEQVS6dWvavN57L1wi2Xbb0EAu7XvmJ08Ot/f16wfDh6tvAylNSuzJUWJP\nWEMT+xtvhJLk5puHkl/HjgUMroV75JHQIPGyyxrfU99jj4WS+u9/DzvvnGx8TTF/frhl7v77wy1z\nhxwS9qm0bvcTaSgl9uQosScs38T+/fdw002h05lTT4Uddih8bBJ6p7vpJthrr9BivkOH/D63eHHo\nunfKlNB17XrrFTbOxvrhh1A1f/fd4Za44cNDSb59+7QjE6mbEntylNgTlk9i/9//QqmxT5/w8JJV\nVmmm4ASAOXPgxhvhpZfgqKPCU+Hqalk+axaMGhUaNP7hD7DSSs0Xa2NVV4ce7e69F95/P8Teuzf0\n6hX+evcOnfioRC/FQok9OUrsOZjZ+sC5wKruPszMNgZOAVYHnnb363J8LmdiX7wYxowJnZCcfHJo\ncCXpef/90Fhx0aLQYHGjjUKjuEWLav5mzQonAQcfHP5KMREuWRJu/fv44/D30Ufh74cfQoLP/G2w\nQbg+365d2hFLS6TEnhwl9nqY2T3uPiw23AoY6+6H55i+1sT+/vuhB7nu3cO91XV1MSrNxx2efRau\nvz7cNrbSSqF6PvO/Qwf4xS9g663TjjR5c+fWJPmPPoIPPwyd7PTrF+6R33LL0jyRkdKkxJ6csk/s\nZnYzsA/wpbtvHhs/GLgaaA3c6O6X5fj8ssRuZj8HfgOMc/c7c0y/XGL/8Ue4/fbQg9xvfwu77aaD\npRSv774L1+jvvz90ZXvAAWGf1TV6KTQl9uS0hMS+M7AAuDWT2M2sNfAesDswA3gFGAFsC2wDXO7u\nM6NplyuxR+Mecfd9cyxvWWKfOjX0INepE5xxhjqbkdJRXR3u/7/vvvDwob33Dg8f6tQp1DZ16hRq\nNTInqe6han/BgvC3aFF4Ut2KK6a7HlI6lNiTU/aJHcDMegIPxxL79sAodx8cDZ8F4O6jY5/pAlwC\n7AbcBLwIHAC0B95w92tzLMvdnQcfDJ2HHHss7LOPSulSuqZPD7cKTp8equ/nzg0NEJcsCQl+yZKQ\nzFu3DrdrduwYSvgzZsB228Guu4be8lTql7oosScn5SdWp6YbMC02PB3oH5/A3WcDx2d97tl8Zl5Z\nWcnMmaFBUseOFZhVNCVWkVStu27o5jjbDz+EJN+uXUjmbdsu//68eaFjoAcfhMsvhwEDQpJfb73Q\nqU/HjoXvfleKV1VVVaO73pa6tdQS+4HAYHc/Lho+DOjv7iclsKzEn8cuUupmzw5J/t//Dk+tmz8/\n3IHQvn1I8iuvHHoE7NUrVOH36gXrrKMH3LQkKrEnp6WW2GcA8b7HuhNK7YmorKykoqKCioqKpGYp\nUtK6dIGhQ8NfRnV1aKw3f34o3WduyXvsMfjkkzCuR4/weNp+/cL/7FoBKX0quSevpZbY2xAaz+0G\nzAQmAyPc/Z0ElqUSu0gCFi4Mt+G9/nroD3/q1JDct9suJPp11007QkmSSuzJKfvEbmZ3AgOB1YAv\ngQvc/RYzG0LN7W43ufulCS3PR40apRK7SMLmzQst9V95JfQa2LcvnHBCqLKX0pUpsV944YVK7Akp\n+8Te3FRiFym8xYvhnnvC3/77w4gRsMIKaUclTaESe3LUNKUAKisrdc1IpIDatYORI+Gf/wy34R15\nZHigks6pS09VVRWVlZVph1FWVGJPmErsIs3vjTfC0/dWXjk8Urexj+WV9KjEnhyV2AtAJXaR5rXl\nlqG//512gpNOCo/n1fl1aVCJPXlFW2I3sw7A9+6+NO1YGkIldpF0ffopXHxx6ML5jDP0wKVSoRJ7\ncoqmxG5mrc3sUDN71My+JNyO9oWZvWNml5vZBmnHKCLFr0eP8DjeXr1Cl84vvJB2RCLNq2gSO/AM\nsAFwNrCOu6/r7msAOwEvA5eZWa2PSi02qooXSVfbtiGpV1bC3/8Of/4zfP992lFJbVQVn7yiqYo3\ns3buvrieadq6+5LmiqkxVBUvUlwWLoSrroKZM0MVvarmi5Oq4pNTNIm9XCixixQfdxg7Fp58EkaP\nDtX1UlyU2JNTNFXxZtYx+t82el66iEgizMK97kceCaeeGrqpFSlXRZHYzewPwAVm9mdgVeC6lENq\nEl1jFylOe+0F558PF10ETz2VdjQCusZeCEVRFW9mAwkN5JYABwF7uvsx6UbVOKqKFyl+U6fC2WeH\nRH/EEaFEL+lSVXxyiqLEDiwEjnT3pe5+N/B02gGJSPnq2TPcEvfCCzBmTNrRiCSrKErscWa2vbu/\nmHYcjaUSu0jpmDMHTjwx9Du/995pR9OyqcSenGIpscetknYATaVr7CKloXPn0Er+ppvCM9+l+eka\ne/KKscS+l7s/mXYcjaUSu0jpeestOO88uPxy6NMn7WhaJpXYk1OMJXYRkWa12WbhqXDnnguzZqUd\njUjTFGNin5J2ACLS8uyyCxx8MJx5Jsyfn3Y0Io1XdFXxpU5V8SKl7R//gA8/hD/9KfQ5L81DVfHJ\nKboSu5ltZ2b/MrPXzWxK9Pdm2nGJSMvwm99Ahw5wzTVpRyLSOEVXYjez94HTgbeA6sx4d5+aVkwN\noRK7SOlbuBB+/Ws46ijYbbe0o2kZVGJPTpu0A6jFl+7+UNpBNEVlZSUVFRVUVFSkHYqINEKHDnDh\nhXDaadC7d+jQRgqjqqpKtwcnrBhL7LsDhxB6n8s8xtXd/f70osqfSuwi5ePxx+Huu+Haa2HFFdOO\nprypxJ6cYkzstwEbA/9j+ar4o1ILqgGU2EXKy+WXw/ffh/vc1ad84SixJ6cYE/t7wMalmh2V2EXK\nyw8/hG5n990X9t8/7WjKlxJ7coquVTzwAtA37SBERADat4fKyvCwmHfeSTsakfoVY4n9XaA38Anw\nQzTa3X2L9KLKn0rsIuVp0qRwj/v118MqJf9Ei+KjEntyijGx96xtfHPf7mZm6wPnAqu6+7BoXAeg\nCqh090dzfE6JXaRM/eUvsHgxnHFG2pGUHyX25BRdVby7T63tL4U4PnH3Y7NG/wG4u7ljEZHicOyx\n8Mor8PrraUciklvRJfakmdnNZjbLzKZkjR9sZu+a2QdmdmYe89kDeBv4qlCxikhx69ABfvc7uPLK\nUHIXKUZln9iBW4DB8RFm1hr4ezS+LzDCzDYxs8PN7Coz61rLfAYCA4BDgePMdOOLSEu0ww7QqxeM\nG5d2JCK1K5rEbmanmlk/M0u0Nzx3fw6YkzW6H/BhVM2/BLgL2N/dx7n7qe4+08y6mNl1wFZmdqa7\nn+fupwJ3AP/UhXSRluvkk+Hhh+GTT9KOROSniqlL2XWBq4FNomrzSYRb315w99kJL6sbMC02PB3o\nH58gWubx2R9097H1zbyysnLZa3UtK1J+VlsNjjkGrrgC/vY3aFU0RaTSoa5kC6cYW8W3B7YFtgd2\niP7PdfdNmjDPnsDD7r55NHwgMNjdj4uGDwP6u/tJTYtereJFWorq6nC9fdAgGDo07WhKn1rFJ6cY\nzzNXBFYBVo3+ZgIvJbyMGUD32HB3Qqk9EZWVlToTFSlzrVqFh8SMGQNfqUlto1VVVS1XyylNVzQl\ndjO7gdCQbT4wGXgReMnds6+PN2bePVm+xN4GeA/YjXDiMBkY4e5N7ldKJXaRlmXMGPjwQ7joIvUl\n3xQqsSenmErs6wHtgS8IJeoZwNymztTM7iRcq9/QzKaZ2VHu/iNwIvAk4Ra2u5NI6hkqsYu0HCNH\nwrRp8OKLaUdSmlRiT17RlNgBzKwVsCk119c3B74hlNwvSDO2fKnELtLyvPwyXHMN3HQTtCmmJskl\nRCX25BRTiR13r3b3KcDj0d/zwAbAKakG1kAqsYu0LP36wZprhlvgpGFUYk9e0ZTYzewUalrB/0io\nPn8++v+Wuy9NMby8qcQu0jJ9/DGcfjqMHQsrr5x2NKVHJfbkFFOJvScwHhjg7r3c/TB3v9bd3yiV\npJ6hErtIy9OrF+y4I9x2W9qRlBaV2JNXTCX2Vu5e3dRp0qYSu0jLNXs2HH10uN7etbaOqSUnldiT\nU0wl9mfM7CQzWy8+0szamdluZnYrcERKsYmI1KtLFxg2LDyzXSQtxZTYhwDVwJ1m9rmZvWNmnwAf\nAiOAq9z9llQjzJOq4kVaroMOgvfegzffTDuS0qCq+OQVTVV8nJm1A1YHvkuig5rmpKp4EXn6abjn\nnlAlr37k86Oq+OQU5S7n7ovdfWapJXUREYBddw0J/emn045EWqKiTOwiIqWsVSs44QS44Qb4/vu0\no5GWRom9AHSNXUQ22wz69oX77ks7kuKma+zJK8pr7KVM19hFJGP6dDjxxNBpzaqrph1NcdM19uSU\nTIndzP5gZmumHYeISL7WXTc8r33cuLQjkZakZBI7cCUw2Mx6pB2IiEi+Dj8c/t//g5kz045EWopS\nSuznANOATdIOREQkX507w4EHhie/iTSHoknsZnaGmY00s71yTDIecOBnzRhWo6jxnIjEDRsWOqx5\n9920Iyk+ajyXvKJpPGdm44BflnrLMzWeE5HaPPJIuK/9yivB1ETsJ9R4LjlFU2IHnnV3N7NOZjbc\nzPYysy5pByUikoQhQ2DOHJg8Oe1IpNwVU2I3AHefS6h2PwjolmpEIiIJad0afvWr8ICY6qJ+RqWU\numJK7Beb2b/MbBQwFHjX3acAmFn/dEMTEWm67beHlVeGJ59MOxIpZ8V0jf0U4GVgALAd0B9oCzwL\ndHD3A1MML2+6xi4idXn7baishFtvhRVWSDua4qFr7MkpmhK7u//F3V9y96uBi9x9A0IL+LuAFdON\nrmHUKl5EcunbFzbaCB57LO1IioNaxSevaErscWZ2LiGpX+Xuz5nZAHd/Ke248qESu4jU5/334bzz\n4PbboW3btKMpDiqxJ6doSuxx7n4xcCCwsZn9tlSSuohIPjbcENZfX9fapTCKMrGb2YYe3ABMNrPt\n045JRCRJhx0Gd94JS5emHYmUmzZpB5DDWDP7BlgKvBmNezHFeEREErX55rDmmvDMM7DHHmlHI+Wk\nKEvswCHuvi/hXvYHgekpxyMikriRI8N1dt3XLkkqysTu7p9G/5e4+6vufn1zx2Bm65vZjWZ2TzRc\nYWbPmdm1ZjawueMRkfLzs5/BSivBpElpRyLlpCgTezFw90/c/djYqGpgPtAe1SCISALMwrX2224D\n3UwjSSn7xG5mN5vZLDObkjV+sJm9a2YfmNmZeczqOXffGzgLuLAgwYpIizNgQGhApz7kJSlln9iB\nW4DB8RFm1hr4ezS+LzDCzDYxs8PN7Coz65o9k9jN6XMJpXYRkSZr1SqU2seNU6ldklH2id3dnwPm\nZI3uB3zo7lPdfQmhd7v93X2cu5/q7jPNrIuZXQdsZWZnmdkvouFbgb8171qISDkbOBDmzYP//jft\nSKQcFOvtboXWDZgWG55O6Jt+GXefDRyf9bl/5TPzePeIFRUVVFRUNCZGEWkhWrWCQw8N19q33jrt\naJpHVVWVut4ukKLsUjZpZtYTeNjdN4+GDwQGu/tx0fBhQH93PymBZalLWRFpsB9/hMMPh/PPD/3J\ntzTqUjY5ZV8Vn8MMoHtsuDsJtnTXQ2BEpKHatIHhw+GOO9KOpHnpITDJa6kl9jbAe8BuwExgMjDC\n3d9JYFkqsYtIo/zwQ6iSv+KK0Jd8S6ISe3LKvsRuZncCLwAbmtk0MzvK3X8ETgSeBN4G7k4iqWeo\nxC4ijdG+PRx4YMsqtavEnrwWUWJvTiqxi0hTLFwYupq95hro+pMbb8uXSuzJKfsSexpUYheRxurQ\nAfbdF+6+O+1ImodK7MlTiT1hKrGLSFPNmQNHHAFjxkCXLmlH0zxUYk+OSuwiIkWmc2fYfXe49960\nI5FSpMReAKqKF5GmGj4cHn0U5s9PO5LCUlV88lQVnzBVxYtIUkaPhm7dQsc15U5V8clRiV1EpEiN\nGAH33w/ff592JFJKlNgLQFXxIpKEHj1giy1ClXy5UlV88lQVnzBVxYtIkt57L/Qff/vt0LZt2tEU\njqrik6MSu4hIEdtoo1BynzAh7UikVCixF4Cq4kUkSSNHhm5mly5NO5LkqSo+eaqKT5iq4kUkae5w\n8snwi1/AoEFpR1MYqopPjkrsIiJFzgwOOyxcZ6+uTjsaKXZK7CIiJaBfP2jdGl56Ke1IpNgpsYuI\nlACzcK39tttC1bxILkrsBaDGcyJSCDvvHB7r+tpraUeSHDWeS54azyVMjedEpJCeegoefxyuuirt\nSJKlxnPJUYldRKSEDBoEX3wBb72VdiRSrJTYRURKSJs2oQ/5229POxIpVkrsIiIlZvBg+PBD+OCD\ntCORYqTELiJSYtq1g2HDVGqX2imxF4BaxYtIof385/DGG/DZZ2lH0jRqFZ88tYpPmFrFi0hzGTsW\nvvwSzjgj7UiaTq3ik6MSu4hIiRo6FJ57Dr7+Ou1IpJgosYuIlKhVV4U99oD77087EikmSuwiIiVs\n2DB49NHQI50IKLGLiJS0tdcOD4h56KG0I5FiocQuIlLiDjkE7rsPFi9OOxIpBkrsOZjZ+mZ2o5nd\nEw23MrOLzeyvZvbLtOMTEcno3Tv8TZiQdiRSDJTYc3D3T9z92Nio/YFuwGJgejpRiYjUbsQIuPtu\nqK5OOxJJW9kndjO72cxmmdmUrPGDzexdM/vAzM7MY1YbAs+7++nAbwoSrIhII225JXTsCM8/n3Yk\nkrayT+zALcDg+Agzaw38PRrfFxhhZpuY2eFmdpWZda1lPtOBudFrnROLSFExC6X2O+8E9ZHVspV9\nYnf354A5WaP7AR+6+1R3XwLcBezv7uPc/VR3n2lmXczsOmDrqER/P7CXmf0VqGrOdRARyceOO8L8\n+fDmm2lHImlqk3YAKekGTIsNTwf6xydw99nA8VmfO5Y8xPs9rqiooKKiojExiog0SKtWMHx4KLVv\nuWXa0dStqqpKz9QokJaa2AteUaWELiJp2GsvGDMGpk6Fnj1TDqYOmWOkEnzyWsRDYMysJ/Cwu28e\nDQ8AKt19cDR8NlDt7pclsCw9BEZEUnXDDbB0KRyfXedYxPQQmOSU/TX2HF4F+phZTzNrBwwHEuu3\nSY9tFZE0DR4c7mn/8ce0I6mfHtuavLIvsZvZncBAYDXgS+ACd7/FzIYAVwOtgZvc/dKElqcSu4ik\n7sQT4dBDYYcd0o4kPyqxJ6fsr7G7+4gc4x8HHi/EMisrK3WNXURSNXgwPPFE8Sd2XWNPXtmX2Jub\nSuwiUgwWLgx9yI8bB506pR1N/VRiT05LvcYuIlLWOnSA7beHp59OOxJpbkrsBaDGcyJSDIYMCdXx\nxUyN55KnqviEqSpeRIpFdTWMHAl//CP06ZN2NHVTVXxyVGIXESlTrVqFDmuKvdQuyVJiLwBVxYtI\nsdhrr3CdfcmStCOpnarik6eq+ISpKl5Eis2pp8LQoTBwYNqR5Kaq+OSoxC4iUuZKoRGdJEeJvQBU\nFS8ixWSXXeB//4Nvvkk7kp9SVXzyVBWfMFXFi0gxuvxy6N49dFpTjFQVnxyV2EVEWoAhQ+Dxx0Hl\njvKnxC4i0gJsuim0aweff552JFJoqopPmKriRaRYVVeHe9uLkarik1OkX3FpU+M5ESlGxZjU1Xgu\neSqxJ0wldhGRhlOJPTlFeP4mIiIijaXELiIiUkaU2EVERMqIEruIiEgZUWIXEREpI0rsBaDb3URE\n8qPb3ZKn290SptvdREQaTre7JUcldhERkTKixC4iIlJGlNhFRETKiBK7iIhIGWmTdgDFyszWB84F\nVnX3YWa2EzCSsM36uvuOqQYoIiJSC5XYc3D3T9z92NjwJHf/DfAIMCa1wEqIbvmroW1RQ9uihraF\nFELZJ3Yzu9nMZpnZlKzxg83sXTP7wMzObMAsDwXuSDbK8qSDVg1tixraFjW0LaQQyj6xA7cAg+Mj\nzKw18PdofF9ghJltYmaHm9lVZta1thmZ2XrAt+6+sNBBZzT0h1/f9HW9X9t72eMaOpwkbYvGz1vb\nIv/ptS3ye7/Yt0VLVvaJ3d2fA+Zkje4HfOjuU919CXAXsL+7j3P3U919ppl1MbPrgK1jJfqjgZub\nL3r9UJsyb22L/KfXtsjvfW2Lut8rpm3RkrWInufMrCfwsLtvHg0fBOzl7sdFw4cB/d39pASWVf4b\nVESkANTzXDJaaqv4giVf7ZgiIpKmsq+Kz2EG0D023B2YnlIsIiIiiWmpif1VoI+Z9TSzdsBw4KGU\nYxIREWmysk/sZnYn8AKwoZlNM7Oj3P1H4ETgSeBt4G53f6eB863zdjkz29jMXjSz783stCTWpVTk\nusWwJaht3aOGmBPM7H0ze8rMOqUZYyE1dP3N7OzoN/Sume2ZTtTJSWr9zexnZjYleu8vzb0eUtrK\nPrG7+wh37+ru7d29u7vfEo1/3N03cvcN3P3Shswz1+1yWZN9A5wEXJHAapSan9xi2ILUtu5nARPc\nfUPg6Wi4XOW9/mbWl1Bb1jf6zDVmVurHpKauf6aNzrXAMe7eh1C72OJ+T2a2kpmdnXYcpajUf0Rp\nqfV2ufgE7v6Vu78KLEkjwDTluMWwRcix7vsBY6PXY4GhzRpUM2rg+u8P3OnuS9x9KvAh4bdVshJY\n//5mtg6wsrtPjqa7lTLeZ+qwNdAu7SBKkRJ743QDpsWGp0fjRGqzlrvPil7PAtZKM5gU5Fr/rizf\naLVcf0cNXf/s8TMoz+2Sk5n1Bj5IO45SpcTeOLpXXRrFQ8cRLXb/yWP9y3rbtPTvH8DMWpnZlVnj\nzjez/czsnGhUd2B9oGc5t0kpFCX2xtHtctIQs8xsbYComvXLlONpbrnWP/t3tG40rtw0ZP2nR+PX\nzRpfFtvFzDoDvwMGxsbtTugs7SGgrZnt7O5VwDxgDUB9gzSQEnvjNOR2Oe2U8hBwRPT6COCBFGNJ\nQ671fwg4xMzaRY9J7gNMruXzpa5B6+/uXwDzzKx/1JjucMpkn3H3Oe5+JSFpZ+wAvBa9fh0YFE37\njrvv6+4tsr1OU7TUnueaxN1/NLPM7XKtgZvc/R0z+3X0/vXRGforwCpAtZmdQniO+4LUAm8m0S2G\nA4HVzGwacEHmboRyF1v31TPrDowGxpvZMcBU4OD0Iiyshqy/u79tZuMJt5z+CJzgJd7HdYLrfwLh\n8dArAo+5+xPNuR7NbE1gUfR6IbB2irGUhRbRV7yIiBQPM5vo7rtGr68h9CXybHQv/37ufmK6EZY2\nVcWLiEiaZgEdoterAF+lGEtZUGIXEZE0TQK2iF5vB7yYYixlQYldRESahZl1MLNTgU3M7Hdm1gF4\nBlgjepy2u/tT6UZZ+nSNXUREpIyoxC4iIlJGlNhFRETKiBK7iIhIGVFiFxERKSNK7CIiImVEiV1E\nRKSMKLGLiIiUESV2ERGRMqLELiIiUkaU2EWawMyuih7Jmxl+0sxuiA3/OepCM8llJvroXzNb1cx+\nExvuaWZT8vzsimZWFT03vCkxtDOzZ82sdVPmIyJK7CJNNQnYAcDMWgGrAX1j728PPJ/wMpPuB7oz\n4fnfjXE0cF9Tn6Pu7ouBp4HhTZmPiCixizTVi4TkDbAp8BYw38w6mVl7YBPgdTN7wMxeNbO3zOy4\nzIfN7FIzOyE2XGlmp5nZYWb2spm9bmbXRScNy8k1TVTifsfM/hkt70kzWyF673wze9fMnjOzO8zs\nNOBSoHc0n8sIJw6ta/t8LQ4FHqxrA0XxvGtmt5jZe2Z2m5ntbmaTzOx9M9sumvQBYGR9G1xE6qbE\nLtIE7j4T+NHMuhMS/IvA5Oj1tsAUd18CHOXu2xIeS3mymXWJZnE3cHBslsOAl6NxO7j71kA1WQnP\nzDauZ5oNgL+7+2bAXODAKIEeQHhE5pAoPgfOAj5y963d/UzAgD7Zn89edzNrB/Ry989i4/Y2s5HR\n64vNrFv0Vm/gCmDj6O8Qd98JOB04J5rmf9H2EZEmaJN2ACJl4AVCdfwOwJVAt+j1t4SqeoBTzGxo\n9HpdQuJ82d3/a2Zrmtk6wJrAHGBL4GfAq9Gl6xWBL7KWuVs903zi7m9Gr/8D9ARWBx6Iqr0Xm9nD\nhCRem9o+n211QtKPGwRk2hhs6e4zzKxnNL//AZjZ/wjV7hBqOHoCuPtSM1tsZh3cfWGOuESkHkrs\nIk33PLAjsDkwBZhGKIl+C9xsZhWERDzA3b83s4lA+9jn7wEOAtYmlOANGOvu55BbfdP8EHu9lJD4\nM5+jltf5fj7uOyC7in5Td38vugwRn0f8dTWwOPY6fhxqD3xfR1wiUg9VxYs03QvAvsA3HswBOhGq\n418AVgHmREl9Y2BA1ufvBkYQkvt4Qmn2IDNbA8DMupjZelmfyWeabM8DPzez9mbWEdiHUBU/H1i5\noSsdrWfrqEoeM1spNp/+wH/NbGC+8zOz1YCv3X1pQ2MRkRpK7CJN9xahNfxLsXFvAnPdfTbwBNDG\nzN4mNFR7Mf5hd38b6AhMd/dZ7v4OcB7wlJm9ATxFKM3HPlL/NFkxuru/CjwUxfYYoXbh2yjG581s\nSqzx3E8+n2PdnwJ2jl73B1Y1s32A7kAHQum7vvllXu8KPJJjOSKSJ2viXSoiUkIy16+j0vWzwHHu\n/t8mzG9r4FR3/6WZnQc85+7PNnJe9wFnuvuHjY1HRFRiF2lp/mlmrxMaxN3blKQO4O6vA89Et9r1\nJqs2Il9m1pbQsE9JXaSJVGIXEREpIyqxi4iIlBEldhERkTKixC4iIlJGlNhFRETKiBK7iIhIGVFi\nFxERKSNK7CIiImVEiV1ERKSM/H8QwCSc4ZrKywAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106420940>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.sed.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$[0.11085907,~0.13624282,~0.16743878,~0.20577778,~0.25289539,~0.31080167,~0.38196693,~0.46942713,~0.57691338,~0.70901111,~0.87135567,~1.0708728,~1.316074,~1.6174197,~1.9877656,~2.4429107,~3.0022718,~3.6897117,~4.534557,~5.5728492,~6.8488825,~8.4170933,~10.344382,~12.712969,~15.623899,~19.201353,~23.597949,~29.001248,~35.641758,~43.802766,~53.832428,~66.158613,~81.307164,~99.924329,~122.80433,~150.92323,~185.48061,~227.95071,~280.14532,~344.2911,~423.12454,~520.00874,~639.07682,~785.40831,~965.24579,~1186.2612,~1457.8833,~1791.6996,~2201.9508,~2706.1386] \\; \\mathrm{\\mu m}$"
      ],
      "text/plain": [
       "<Quantity [  1.10859065e-01,  1.36242824e-01,  1.67438785e-01,\n",
       "             2.05777785e-01,  2.52895390e-01,  3.10801666e-01,\n",
       "             3.81966929e-01,  4.69427133e-01,  5.76913383e-01,\n",
       "             7.09011107e-01,  8.71355674e-01,  1.07087280e+00,\n",
       "             1.31607401e+00,  1.61741974e+00,  1.98776557e+00,\n",
       "             2.44291068e+00,  3.00227185e+00,  3.68971175e+00,\n",
       "             4.53455699e+00,  5.57284917e+00,  6.84888248e+00,\n",
       "             8.41709327e+00,  1.03443824e+01,  1.27129691e+01,\n",
       "             1.56238987e+01,  1.92013534e+01,  2.35979494e+01,\n",
       "             2.90012482e+01,  3.56417578e+01,  4.38027665e+01,\n",
       "             5.38324276e+01,  6.61586128e+01,  8.13071645e+01,\n",
       "             9.99243290e+01,  1.22804326e+02,  1.50923230e+02,\n",
       "             1.85480610e+02,  2.27950706e+02,  2.80145316e+02,\n",
       "             3.44291096e+02,  4.23124542e+02,  5.20008738e+02,\n",
       "             6.39076821e+02,  7.85408309e+02,  9.65245791e+02,\n",
       "             1.18626124e+03,  1.45788332e+03,  1.79169958e+03,\n",
       "             2.20195083e+03,  2.70613864e+03] micron>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.sed.wavelength"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 50)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.sed.flux.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with Model Results: Images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each ModelResults instance will have an attribute ``.images`` which is a collection containing the different image outputs. This looks like a Python ``dict`` but with some convenience customizations for use with MCFOST. For instance you can index it with either strings or floats for the wavelength interchangably. \n",
    "\n",
    "Each item in this dict is an ``astropy.io.fits`` HDU containing the image data and its header. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$[0.8,~1] \\; \\mathrm{\\mu m}$"
      ],
      "text/plain": [
       "<Quantity [ 0.8, 1. ] micron>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.images.wavelengths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_values(['0.8', '1.0'])\n"
     ]
    }
   ],
   "source": [
    "print(model.images.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<astropy.io.fits.hdu.image.PrimaryHDU object at 0x10e524b70>\n",
      "<astropy.io.fits.hdu.image.PrimaryHDU object at 0x10e524b70>\n",
      "<astropy.io.fits.hdu.image.PrimaryHDU object at 0x10e524b70>\n"
     ]
    }
   ],
   "source": [
    "print(model.images[0.8])\n",
    "print(model.images['0.8'])  # This is the exact same \n",
    "print(model.images[' 000.8000'])  # So is this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SIMPLE  =                    T / file does conform to FITS standard             \n",
       "BITPIX  =                  -32 / number of bits per data pixel                  \n",
       "NAXIS   =                    5 / number of data axes                            \n",
       "NAXIS1  =                  101 / length of data axis 1                          \n",
       "NAXIS2  =                  101 / length of data axis 2                          \n",
       "NAXIS3  =                    1 / length of data axis 3                          \n",
       "NAXIS4  =                    1 / length of data axis 4                          \n",
       "NAXIS5  =                    4 / length of data axis 5                          \n",
       "EXTEND  =                    T / FITS dataset may contain extensions            \n",
       "COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy\n",
       "COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H \n",
       "WAVE    =                   1. / wavelength [microns]                           \n",
       "CTYPE1  = 'RA---TAN'                                                            \n",
       "CRVAL1  =                   0. / RAD                                            \n",
       "CRPIX1  =                   51                                                  \n",
       "CDELT1  =            -5.89E-05 / pixel scale x [deg]                            \n",
       "CTYPE2  = 'DEC--TAN'                                                            \n",
       "CRVAL2  =                   0. / DEC                                            \n",
       "CRPIX2  =                   51                                                  \n",
       "CDELT2  =             5.89E-05 / pixel scale y [deg]                            \n",
       "BUNIT   = 'W.m-2.pixel-1'                                                       \n",
       "FLUX_1  = 'I = total flux'                                                      \n",
       "FLUX_2  = 'Q       '                                                            \n",
       "FLUX_3  = 'U       '                                                            \n",
       "FLUX_4  = 'V       '                                                            "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.images[1].header"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For now there is a \"plot_image\" function to display images. Maybe eventually this should become a method on a ModelImage object once such a thing would exist? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cpinte/Software/miniconda3/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  if self._edgecolors == str('face'):\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAEbCAYAAAAPs+VNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXFWd9/HPl7AvCQSEkBAkKGEIIJsmjOADI4KRLa4s\nOsjqM4iMjKIsMjqO24CoMMDDqMMqalgUFJBFRFEWAdkEDIFEE00ISUSIhCWQkN/zxz03uanqrq6q\nrr5V1f19v171qrrn3jrnVCf9q9PnnkURgZmZtddq7a6AmZk5GJuZdQQHYzOzDuBgbGbWARyMzcw6\ngIOxmVkHcDA2M+sADsZmZh3AwdhaStK2kh6R9IKkE9tdH7Nu4WDcRSTNlrRPu+vRh1OA2yNieERc\n0N/MJI2UdJ2kF9PnP7zGtVtJuknSc5KekXS+pGH9rYNZGRyMu0ukRyd7IzCt0TdJWr2XU/8PWAJs\nCnwE+B9JE3q59kJgATAK2BnYCzih0bqYtYODcZdKrcTPSPp9ajVeJGkzSTenLoLbJG1YuP40STPT\nuT9Iem/h3K6SHk7nrpZ0laQvF86PlvRjSQsl/UnSv/ZSp18CewMXpLzeLGk7SXdIel7S45IOqvgM\np0h6FFgsabWK/NYD3g98PiJejoi7gZ8CR/TyY9kKuCoiXouIBcAtwPYN/WBXrVvdP1+z/nIw7m7v\nB94FjAcOAm4CTgPeQPZv+8nCtTOBPSNiOPCfwPdTcFkTuA64BNgImAq8l9QCTwHyBuBhYDSwD/Bv\nkvarrExEvBO4E/hEKufP6b23pDr9K/ADSdsU3nYY8B5gw4hYXpHleGBZRMwspP2e3gPsucBhktaR\nNCble3Mv1/YlaOzna9YvDsbdK4DzI+KvETGPLAjeGxG/j4hXyQLsLisujvhRRMxPr68GZgCTgN2B\nYRFxfkS8HhHXAfcXynkbsElEfCUilkXELOAisiDal92B9SLizPTeXwE3Ah8ufIbzIuLpVOdK6wMv\nVKQtBjbopbw7gR3Se+YAv4uIn/ZWOUnDJP2mcPxtSdsWLqn752vWXw7G3W1B4fUrFcdLyIIZAJI+\nmroinpf0PFnQ2gTYHHi6It85hddvBEbn70vvPZ2sD7cvoyvygqy1PLqXsiq9CAyvSBtBFpBXkVrw\ntwA/BtYl+2wjJZ1VI//dgT+m9wv4PxHxZOF83T9fs/5yMB5c1GOi9Ebgu8AngJERsRHweDr9DDCm\n4i1bFl7PAWZFxEaFx/CIOLCO+swDxqZAl3sjqwb/WjcknwJWl/TmQtpOhboXjQTGAhdExNKIeA64\nDNi/Rv6TgdvS6x2BJ2pcC738fM1awcF4aFiPLOg9C6wm6WiyljHAb4HXJZ0oaXVJU8i6JnL3k91c\nOyX1xQ6TtIOkt9YoLw9a9wEvA6dIWkPS3sCBwJX1VDoiXgKuBb4kaV1Je5L13V7Rw7XPArOAj6c6\nbggcSdbH3Jv9CucPAH4p6eB66mbWag7Gg0tUvA6AiJgGfJMs8M4nC8R3pXNLyW5UHQs8TzZ87Ebg\ntXT+dbIAujPwJ+CvZK3syu6DqnpExGtkwfM96X0XAEdExFMNfKYTgHWAhcD3geMjYkULNo0rPi0d\nvr9Q1gzgVeBTPWUqaROy0RcHSzqArNthk/Rc83MVXnf6MEPrIvK2S1ZJ0n3AhRFxebvrMlAkfRjY\nPiLOaHddzMAtYwMk/R9Jo1I3xZFkLedb2l2vAbY7WReIDTGSxqVx49ek4ymSvivpSkn7tq1ebhmb\npI8BXybrW/4jcHpENDs+16wrSLomIj5UON4Q+EZEHNeO+rhlbETE/0bEqIjYICJ2diC2biDpEkkL\nJD1WkT5Z0nRJMySd2kCW/052X6MtHIzNrFtdSjY8cYW0MNQFKX0CcHiakn+EpHMkja7MRJmzgJsj\n4pEyKt6T3hZn6WiS3Ldi1kYR0fSY62Z+f3sqLyLulLRVRfJEYGZEzE5lXQlMiYgzSUMiJY0Evgbs\nnEbivEQ2zX+4pDdHxHcarV8rdGUwNrPutnx55TIkK91xxx3ccccdK46/9KUvNZL1GFad1TmXbNr/\nCmlC0PEV7zu/kUIGQlfewHPL2Ky9+tsyfv311+u+ftiwYb2Wl1rGN0TEjun4A8DkiPhYOv5nYFJE\n9LjSYCdxy9jMSjeAjcCnyabF58aStY47noOxmZWuVjdFPz0AbJNazPOAQ4Fed4fpJB5NYWali4i6\nH72RNBW4BxgvaY6koyNiGXAicCvZjjNXFafPdzL3GZtZw/rbZ7xkSa0lQFa19tpr96u8buFuCjMr\nXTc2Ageag7GZlc7BuJqDsZmVzsG4moOxmZXOwbiag7GZlW4Ah7Z1LQdjMyudW8bVHIzNrHQOxtUc\njM2sdA7G1RyMzax0DsbVHIzNrHQOxtUcjM2sdA7G1dq+UJCkYZIelnRDOh4p6TZJT0n6edok0MwG\nkeXLl9f9GCraHoyBk8hWV8q/Kk8DbouI8cDt6djMBpFWrNo22LQ1GEvaAtgfuAjIV2U6GLg8vb4c\neG8bqmZmA8jBuFq7+4zPAT4LDC+kbRYRC9LrBcBmpdfKzAbUUAqy9Wpby1jSgcDCiHiYla3iVUT2\nL+Z/NbNBxi3jau1sGb8dOFjS/sDaZNtkXwEskDQqIuZL2hxY2MY6mtkAGEpBtl5taxlHxOciYmxE\njAMOA34ZEUcA1wNHpsuOBH7Srjqa2cBoZ8tY0jhJF0m6Jh3/g6T/kXSNpONbXmCdOmE0RS7/qZ8J\n7CvpKeCd6djMBpF2Dm2LiFkRcVzheHpEfJxs89I9Wl5gnToiGEfEryPi4PT6uYh4V0SMj4j9ImJR\nu+tnZq3Vog1JL5G0QNJjFemTJU2XNEPSqfXUR9JBwI3ATf36YP3QEcHYzIaWFnVTXApMLiZIGgZc\nkNInAIdL2k7SEZLOkTS6l/rcEBH7Ax9pzSdsXLuHtpnZENSKvuCIuFPSVhXJE4GZETEbQNKVwJSI\nOBO4IqWNBL4G7CzpNOC3wPuBtYCf9btiTXIwNrPS1QrG9957L/fee2+zWY8B5hSO5wKTKsp+Dqi8\nUffrZgtsFQdjMytdrWA8adIkJk1aGT/PO++8hrJuvlbt5WBsZqUbwAWAngbGFo7HkrWOO56DsZmV\nbgAnfTwAbJP6kueRDVc7fKAKayWPpjCz0rVoaNtU4B5gvKQ5ko6OiGXAicCtZKtBXhURT5TyofrJ\nLWMzK12LRlP02OKNiJuBm/tdQMkcjM2sdF6bopqDsZmVzsG4moOxmZXOwbiag7GZlW4o7W1XLwdj\nMyudW8bVHIzNrHQOxtUcjM2sdA7G1RyMzax0DsbVHIzNrHQOxtUcjM2sdA7G1RyMzax0HtpWzcHY\nzErnlnE1B2MzK52DcTUHYzMrnYNxNQdjMyudg3E1B2MzK52DcTUHYzMrnYNxNW+7ZGalW758ed2P\nVpM0TtJFkq4ppK0n6XeSDmh5gXVyMDaz0rViD7x+lD0rIo6rSD4FuKrlhTXAwdjMSteiDUkvkbRA\n0mMV6ZMlTZc0Q9KpfdVF0r5km5f+td8frB/cZ2xmpWtRi/dS4Hzge3mCpGHABcC7gKeB30m6Hngr\nsCtwdkTMq8hnL2A9YALwiqSbog2d2g7GZla6Fu0OfaekrSqSJwIzI2I2gKQrgSkRcSZwRUobCXwN\n2FnSqRHx7yn9SOCv7QjE4GBsZm1QK9498sgjPPLII81mPQaYUzieC0yqKPs54Pge6nR5PQWkYN6X\n5RGxqJ78cg7GZla6WsF4p512Yqeddlpx/L3vfa/Xa3vKuvla1e0ZoLKro9LqwNhGMnUwNrPSDeCq\nbU+zahAcS9Y6bqUnImLnWhdIarhp72BsZqUbwG7ZB4BtUl/yPOBQ4PAWl7F7i65ZhYe2mVnpWjS0\nbSpwDzBe0hxJR0fEMuBE4Fay4WpXRcQTLa7+v0t6n6RRFfXZK+9PjogljWbqlrGZla5Foyl6bPFG\nxM3Azf0uoHfrACOBr0gaAywC7iVrlR8LnN1Mpm0LxpLGko0P3JSs0/27EXFe+ma5CngjMBs4pNG7\nkmbW2bp5bYqIODm9vBhA0oZk3RL7AH9qNt92dlMsBT4VEduTfZBPSNoOOA24LSLGA7enYzMbRNo5\nHbqVUlfFJhFxC1nj8ofN5tW2YBwR8yPikfT6ReAJsjGCBwP5eL/Lgfe2p4ZmNlDauVBQi30A2FLS\nPwF/ox83CzviBl6687kLcB+wWUQsSKcWAJu1qVpmNkAGS8sYWDMifgmsFxGLgb83m1Hbb+BJWh/4\nMXBSRCyWtOJcRISkjv/XMLPGdEGQrdd0SXcCMyStDrwFuLGZjNoajCWtQRaIr4iIn6TkBZJGRcR8\nSZsDC9tXQzMbCIMlGEfEzZKmkXWnbgr8d7N5tXM0hcjuRk6LiHMLp64HjgTOSs8/6eHtZtbFBksw\nlrQ18ExENB2Ec+3sM94D+GfgnyQ9nB6TgTOBfSU9BbwzHZvZIDKI+oxPJi1EJOkdkt7RbEZtaxlH\nxF30/mXwrjLrYmbl6oIgW6/7gXGS/pyW9Hxfsxl1xGgKMxtaBtHQtrHAq8CnJf0K2K3ZjNo+msLM\nhp5B1DL+E/DjiPihpE2A9zebkVvGZla6QdRnfBWwQ3o9jn7Mi3AwNrPSdXMwTqvF5d4H/EOaL7EG\ncHez+bqbwsxK14lBtgEfLbweTTYN+hKyBc8WAr9sJlMHYzMrXTcH44hYWji8ERgVEYdI2oB+xNRe\n3yjp5N7OFbwYEd9ptnAzG5q6ORhX2Ac4XNKawPci4rvNZlSrz/gzwPo1Hhuka8zMGtLOoW2Sxkm6\nSNI16XhvSXdK+h9JezWY3d8i4p1kq02+Jun0ZutVq0n9/Yj4z1pvlrReswWb2dDVzpZxRMwCjsuD\nMbAcWAysReObl64jadeIeAi4TNLBzdar12AcEZ/t6831XGNmVqkVwVjSJcABwMKI2LGQPhk4FxgG\nXBQRZ/WR1Z0R8RtJmwLfIlumoV47ArtI+mo6fk3SYmCLiLiigXz6Htom6d8kjVDm4rSGxLsbKcTM\nrKhFQ9suBSYXEyQNAy5I6RPI+nO3k3SEpHMkje6hLnkhi8hax424Hrg2It5DtnLb14G3Ap9uMJ+6\n7vwdExHnpgA8EjgCuIJs91Uzs4a1aEPSO9PGFEUTgZkRMRtA0pXAlIg4kyxukfbZ/Bqws6TTgCeB\ndwMbAufXKlPStsDyiJiR6nBPoT6vko0zvlvStY1+nnqCcb7a+wFk6w4/XlwA3sysUbWC8fTp05k+\nfXqzWY8B5hSO55JWVSuU/RxwfMX7rqsz/z8Ce0vaj6yv+XcR8UDlRRHxx7prnNQTjB+U9HNga+B0\nScNTJczMmlIrGG+77bZsu+22K45/+tOfNpR187WqI/OIZcAv0gNJEyV9nKzL90ngjnRNw+rqpiDb\nn+6PEfGSpI2Bo5spzMwMGMjV2J4mW0ktN5bGR0jULSLuJ1tGM+/CODaNOX4auDUiXqo3r1qTPnZj\n5bdMAFun7gkxwN8+Zja4DeDQtgeAbVJf8jzgUPqxY3MjIuJJstYx6UbhgWQLCdWlVsv4m9QOuv9U\nbyE2NGywwQYrXm+00UYArLvuugCstVZ2k7piw1kAXn31VQBefvllAJ5//vkV1yxevHgAa2zt0qKh\nbVOBvYCNJc0BvhARl0o6kWyAwTDg4oh4ot+FNSgi5tFAIIba44z37m+FzMx60qLRFD22eCPiZuDm\nfhdQsj77jNMsu08DW0bExyRtA2wbEU1tR22Dx9ixWdfc6NHZ0M3hw4evOLfmmmsCMGzYMABWW633\nIe15/+Hrr78OwGuvvbbi3AsvvADAvHnzAJgzZw7W/QbR2hQtU896xpcCrwFvT8fzgK/2frmZWW1d\nvp7x+ul5jTTJpCXqGU3xprQ83GEAaURFq8q3LrLeetlSJDvskG1s8IY3vAFY2R9cbP3m/0cqn3uS\n/8Llz3k/M6xsbedl5a3wxx9/HICXXqr7ZrV1kE4MsvWQdAqwSQrC/5UeH2tF3vUE41clrVOozJvI\nNuAzM2tKF2w02pv70mMp8EFauFtSPcH4i8AtwBaSfgjsARzVqgpYZyuOkJg4cSKwcqTE6qtn/33y\nVm9PLePejosqW0nF4/yXdo011gBg7bXXBmD99dcH4P77719xrUdedI9ubRkDLwFHRcS3gauK3RSS\ndgYejYimvmn6DMYR8XNJDwG7p6STIuKvzRRmZgbdG4zT1Ofi9OfVJZ0D/A74DVlD9ZJm8q5n1bbv\nA8si4sY0gmI9SU3t8WRmBt19A68HXyZb8e00YONmM6mnm+JO4D5JnyZbhOMzQD1bMtkgkHdNAGyy\nySZA790Txa6Iem7cVaq8kVd8nXdX5EPl8roU63f77bfXXZa1V5cE2Xo8C7wWETcBN/Uno3q6Kb4j\naRrZjqfPArtGxDP9KdTMhrZBFIwnA5+V9BzZGhW/SutVNKyeSR9HAF8g2576LcBNko6OiEeaKdC6\nw/bbbw+sHFIGK2+i5S3hyokcxePerulpOnQub/0W77RXtojz4zzfYv3yOv/hD3/o6+NZm3XxaIpK\nd0TEJyWtS7ao/FtJCwc1qp5uig8Ae0TEQmCqpOuAy4CdmynQzGwQtYxD0tsiIr+B95tmM6qnm+K9\nFcf3S5rY2/U2OIwbNw5YOaEDqlu7vT0XX9ea/FHZR1zZP1x83dtzscy8zm4Zd75BFIz3ApD0H8Ar\nwK8j4oJmMqq1hOapEXGWpJ62IQngk80UaGY2iILxj8i20bs7dVVMaDajWi3jaen5wR7ODZqfpPUs\nn1RRbHnmoygqW8J5f25+HqpHWtSzUFD+XPxFXbYs2zQhX0So8tr8fLHO1vm6NRj3sAfeXfm5iHiZ\nVccgN6TWEpo3pOfLUiVGpEp4mpOZ9Uu3BmPq3AOvGfWMpngb2YyS4el4EXBsqypgnSlv2eYjKKC6\nJVz5XGwZ568rrym2kCuXzsyfi63dyhZ25bXFPuglS5Y0/kGtLbo1GLd7D7xLgBMi4s5U+J4p7S3N\nFGhmNliGtrVyD7x6VhxalgfiVPhdQFORv16SJkuaLmmGpFMHsiwzK187p0NLGifpIknXpOPVJH1V\n0nmSPlrH+7frKT0inoyI70TE+cC9ZHvg1a2elvGvJX0HmJqOD01pu6YKPNRIgX1JqyBdALyL7Nvl\nd5Kub8c+VkNZZZdE8XVlt0T+XOzSyHf6yJ/XWWedqmuWLl0KwCuvvAKs3OGjuNNHXo+866JyiFze\nXVG81jpfO7spImIWcFwejIEpZEs9PEt9O0lfJOngiPhbjTJatwdewc5koyf+o4d0aP3GpBOBmREx\nG0DSlWQ/LAdjs0GiRRuSXgIcACyMiB0L6ZOBc8k2JL0oIs7qI6vxwN0R8b8pQPe1ENp/A+MlbZze\n93wf19elnkkfe7eioAaMAYobnc0FJpVchyHvxRdfBFbd167yJlxlizhvBcPKdYfz3UHyhX0OOeSQ\nFddcffXVADz77LNAfbt2VPY1FlvGeZ2t87WoZXwpcD7wvTyht7+syaYp7wqcnVqtRXPJtpaDbIRE\nTRFxdaG8t0saCdwVEYv68Vl67zOW1Gd/Rz3XNKE7b7OaWd1a0Wec7mVVtkpX/GUdEUuBK4EpEXFF\nRHwqIuZJGinp28Au6Z7UtcC7JZ0H3NFX3SUdWjh8BHgYOFTSpyQN7+VtfarVMv6GpHmA6DlAimz/\np1bvEv00MLZwPJb6+nGshf72t6w7LG/Rwsqp0X1N/ii+ruzzve6661Zck6f1NmSupzLysvP+5mL/\ncl5n63y1guzs2bOZPXt2s1n3+Zd1RDwHHF/xvuMaKOPiNDP5VeBFsrWMn0/PG1LdpVuXWsF4PvDN\nPt7/VDOF9uEBYBtJW5HtRH0ocPgAlGNmbVJraNuWW27JlltuueL417/+dSNZl/GX9THAbcD+wLMR\ncWsrMq01A2/vVhTQqIhYJulE4FayDviLPZKifPnoh4ULF65Iy/uEiyMiYOXIhuIvWN6Xm7dcX3jh\nBWDVfuHKayqnPBfzrvTqq69W1S+vs3W+ARxNUcZf1jdFxIvADyRtLul44C9pgfmm1TOaonQRcTNw\nc7vrYWYDYwCDcRl/WV+WbgrmLYVXgLen/uez0/Z0DevIYGxmg1uLhrZNJVvCcmNJc4AvRMSlJfxl\nvTPwOiv7iZ8n66e+kCwwN6WetSnWjoglfaXZ4LLuuusCq3YrzJ8/H4BRo0YBsMEGGwAruxeKky6K\n60sUj3taz7hyTYqeujvy53z4Wl6XfMIIrBxGZ52vFcE4Inps8Zbwl/X7I+LRVmdaz5Sle+pMMzOr\nSzfvDj0QgRhqLy6/OTAaWDdNfc6HuA0H1h2IyljnyG+qbbzxyp3HFyxYAKwcVpaf23DDDYGVrWmo\nXpGtpx2kK3f2yFvGxeFqL7/8MgCLFmXj6fPha3kLebPNNltxbX6T0DpfJwbZekl6KCJ27e81lWp1\nU+wHHEU2bq84xG0x8LlGCjEzK+ryVdu2k/RYH9eMaDTTWkPbLgcul/TBiPhRoxlbd5s1axYAu+yy\ny4q0vBWat07nzs1GDOWt1mLLON91I58WXTlRBKpbxPl6xMVpzXnLOE/LW9Z5XYr55XW2ztfNLWOg\nx1XbKjS8smU9oynuknQxMCYiJkuaAPxjRFzcaGFmZtDdwThfxKzV6gnGl5EtyHFGOp4BXA04GA9i\n06dPB1hlJtSmm24KrGyV5iMtFi/OduIq9tlWLq+ZP/fUZ5y3jCufYWXLNx+5kY+YyPMpTvrI62yd\nr5uD8UCpZzTFJhFxFdm4OtLiGwO6uLyZDW7dPJpioNTTMn4xrdsJgKTdgb8PXJWsk/z2t79d8XrP\nPfcEVi4elI+iyFut+RTl4uvKBeOLv1yV+9vlfc75gkTF1/niQfnojHzZzWL9rHsMhiAr6ZPAFaWt\nZwycDNwAbC3pHuANwAdbUbiZDU2DIRgDm5Gtl/wQ2b6gt0Y/Plg9i8s/KGkvYNuU9GTqqjAza0qX\nD20DICLOkPR5Vg4DvkDS1WRTsP/YaH71TIc+BLglIh5PBe8i6Sut3vvOOlN+cw7g7rvvBmDXXbOx\n7PmNvHy1tOJ05N6mJvfUTVFLfn0+7TmfePLQQw9V1c+6xyBpGRMRyyXNBxaQ3VfbCPiRpF9ExGcb\nyaueG3ifj4gXJO0J7EPWHP92o5U2M8sNhht4kk6S9CDwdeBuYIeI+DiwG/D+RvOrp88432TsQOB/\nI+JGSV9utCDrfvnQtTvuuAOA8ePHAzB2bLZ87EYbbbTi2t52BempNVw5+aN4I/D557N7I3PmZJs3\nPPXUQOxnYGXr5CDbgJFkiwb9uZiYWssHNZpZPS3jpyV9l2xd0J9JWrvO95mZ9WgwtIyBdSoDsaSz\nACJiWqOZ1VooaOuI+BPwIeA9ZIsmL0oLCDXUF2KDU95KzZ/zpTVh5SJC+bC3fJeQ4v52lTt95P2/\nxb3s8qUybXDp8CBbr317SNsfOLWZzGp1U1xD1vdxQ0TskydGxDPAM80UZmYG3T2aQtLHgROAN1Us\nGLQBWd9xU2oF42GSzgC2lfRpVm4xAhAR8a1mC7XBqdiKdYvWaunylvEPyRavP5OsFZzHxsUR0fQW\n5bWC8aHA+8i2Ltmg2QLMzCq1MxhLGke21s6IiPhQGin2EbJ4OCEi9qj1/oj4O9ks5MNaWa9awXhy\nRJwpac2I+FIrCzWzoa2dwTgiZgHHSbomHd9FtjrlFOD+vt4v6e6I2EPSi2QbblRkH8ObqVetURHH\npOf3NZOxmVlvWjGaQtIlkhZULvQuabKk6ZJmpB2b6/Vhsi6Ivuq+R3pePyI2qHg0FYihdjCeJmkG\nWZ/xYxWPAdkDysyGhhYNbbsUmFxMkDQMuCClTwAOl7SdpCMknSNpdE8ZSdoS+HtEvNTT+V7e8yFJ\nG6TXn5d0bdqirim1dvo4XNIo4OfAQax6A8/MrGmt6KaIiDslbVWRPBGYmS8AL+lKYEpEnAlckdJG\nAl8jW9rh1Ig4i6wn4JIGq/CFiLimMDv5G2Szkyc283lqzsCLiPmSJgFvSkkzI2JJMwWZmeVqDW2b\nP39+f0bjjAHmFI7nApOKF0TEc8DxFWlfbKKsls5OrjXpYw3gq2TfGH9JyVtKuhT4nFduM7Nm1WoZ\nb7bZZqvs+v3oow31ipZ5ZzCfnbwvcFZ/ZyfXeuPZZHOvx0XErpFtO701sCFZc9zMrCkDOB36aWBs\n4XgsWet4IBwC3ALslxaY34h+zE6u1U1xIDA+Ilb8PZFWbzseeBI4qdlCzWxoG8ChbQ8A26S+5Hlk\n8yUOH6CyXgfWAQ6RlMfSILvP1rBaLePlxUCci4jXge6dy2hmbdeioW1TgXuA8ZLmSDo6IpYBJwK3\nAtOAqyLiiQH6GD8FDgaWAi+mR92jMSrVahk/IenIiLi8mCjpCMDb8JpZ01o0mqLHFm9E3Ew2XXmg\njYmId7cqs1rB+BPAtZKOAR5MabsB6+KJIGbWD12+NkXuHklviYiWzLuoNc54bhrW9k5ge7K+kJ9F\nxO2tKNjMhq5uXrWt4B3A0ZJmAfmOCBERb2kms77GGQdwe3qYmbXEIGkZvyc9By2YFOcdO8ysdINk\np4+/kLWOj0wz/pYDmzabmYOxmZVukATjC4F/JFtgCLLRFBc2m1lbgrGksyU9Ien3aXGNEYVzp6fV\nlqZL2q8d9TOzgTVIgvGkiDgBeAVWTLNeo9nM2tUy/jmwfUTsBDwFnA4gaQLZIO0JZKsuXSjJrXez\nQWaQBOPX0ipxAEh6A/2Yg9GWQBcRtxUmlNwHbJFeTwGmRsTS1AczkyZXQDKzzjVIgvH5wHXAppK+\nRrb/3X81m1nN0RQlOQaYml6PBu4tnJtLtgqTmQ0ig2FoW0R8X9KDZMtnQrZUZ9Oz/QYsGEu6DRjV\nw6nPRcQN6ZozgNciotbq+h391WhmjevwFm9dJJ0VEacCT/SQ1rABC8YRsW+t85KOAvZn5bcKVK+4\ntEVKM7NBZDAEY2A/st2hi/bvIa0u7RpNMZlsqbkpFYvVXw8cJmnNtIPrNtSxQaCZdZdu7jOW9PG0\n717llnQ4ipsoAAAMJ0lEQVSzgaanRrerz/h8YE3gNkkAv42IEyJimqSryVZbWgacEJ34r2Fm/dLl\nv9Y/JFuI6EyyVnA++25xRPyt2UzbEowjYpsa575Gtj+VmQ1S3RyMI+LvwN+Bw1qZbyeMpjCzIaab\ng3EubbP0AWArVsbSiIgvNZOfg7GZlW4wDG0jW1x+EdkSw/3eqNnB2MxKNxhaxrR4cXlPNTaz0rVz\nNIWkcZIuknRNOt5S0nWSLpbUyLC0eyQ1tXZxTxyMzax07QzGETErIo4rJO0A/CgijgV2aSCrdwAP\nSnqqMLyt64a2mdkQ1oogK+kS4ABgYUTsWEifDJwLDAMuioiz+sjqPuCatMXcFQ1UYXKDVa7JwdjM\nSteiFu+lZHMWvpcnpFXULgDeRTZ793eSrgfeCuwKnB0R8yryOQr4j4i4M3VdXFarUEkv0vsyDQEM\nb/iT4GBsZm3QitEUKXhuVZE8EZiZVn1E0pVkM33PJLV6JY0km8uwc+ojvhH4oqQPA7PqKHf9fle+\nBw7GZla6Wi3jRYsWsWjRomazHgPMKRzPBSZVlP0ccHzF+z7UbIGt4mBsZqWrFYxHjBjBiBErNv/h\nL3/5S0NZN1+r9nIwNrPSDeA448qVH8eStY47noOxmZVuAIPxA8A2qS95Htk2bocPVGGt5HHGZla6\nVowzljQVuAcYL2mOpKMjYhlwInAr2eqPV/Vn940yuWVsZqVrRcs4Inps8UbEzWRLXHYVB2MzK90g\nWSiopRyMzax0g2ShoJZyMDaz0jkYV3MwNrPSORhXczA2s9I5GFdzMDaz0jkYV3MwNrPSORhXczA2\ns9J5aFs1B2MzK51bxtUcjM2sdA7G1RyMzax0DsbVHIzNrHQOxtUcjM2sdA7G1RyMzax0DsbVHIzN\nrHQe2lbNwdjMSueWcTUHYzMrnYNxNW+7ZGala8W2S82SNE7SRZKuSccTJF0l6UJJH2h5gXVyMDaz\n0rUzGEfErIg4rpA0GTg/Ik4APtryAuvkYGxmpWvRhqSXSFog6bGK9MmSpkuaIenUOqpzBXCYpK8D\nG/fvkzXPwdjMSteilvGlZK3aFSQNAy5I6ROAwyVtJ+kISedIGt1DXf4aEScCpwPPtuxDNqitwVjS\nyZKWSxpZSDs9faNNl7RfO+tnZgNj+fLldT96ExF3As9XJE8EZkbE7IhYClwJTImIKyLiUxExT9JI\nSd8GdpZ0qqQ3SvoOcDnw9QH6yH1q22gKSWOBfYE/F9ImAIeSfaONAX4haXxEeFCi2SBSq8W7ZMkS\nlixZ0mzWY4A5heO5wKSKsp8Djq943780W2CrtLNl/C3glIq0KcDUiFgaEbOBmWTfdGY2iNTqllhr\nrbUYMWLEikejWQ9EfcvQlmAsaQowNyIerTg1muybLDeX7JvOzAaRARxN8TQwtnA8llVjSscasG4K\nSbcBo3o4dQZZR3mxP1g1surabzoz69kATvp4ANhG0lbAPLJuz8MHqrBWGrBgHBH79pQuaQdgHPB7\nSQBbAA9KmkT1t9oWKc3MBpFWBGNJU4G9gI0lzQG+EBGXSjoRuBUYBlwcEU/0u7ASqN3TEiXNAnaL\niOfSDbwfkvUTjwF+Abw5Kiopya1lszaKiFp/zdYkKTbffPO6r3/mmWf6VV636IS1KVYE1oiYJulq\nYBqwDDihMhCbWffzqm3V2t4yboZbxmbt1d+W8aabblr39QsXLnTL2MxsIHRjI3CgORibWekcjKs5\nGJtZ6RyMqzkYm1npHIyrORibWek8mqKag7GZlc4t42oOxmZWOgfjag7GZlY6B+NqDsZmVjoH42oO\nxmZWOgfjag7GZlY6B+NqDsZmVjoPbavmYGxmpXPLuJqDsZmVzsG4moOxmZXOwbiag7GZlc7BuJqD\nsZmVzsG4moOxmZXOwbiag7GZlc5D26qt1u4KmJlZl25IamY22LhlbGbWARyMzcw6gINxBUknS1ou\naWQh7XRJMyRNl7RfO+uXk3S2pCck/V7StZJGFM51XH1zkianes2QdGq761NJ0lhJv5L0B0mPS/pk\nSh8p6TZJT0n6uaQN213XIknDJD0s6YZ03NH1tWoOxgWSxgL7An8upE0ADgUmAJOBCyV1ws/t58D2\nEbET8BRwOnR0fZE0DLiArF4TgMMlbdfeWlVZCnwqIrYHdgc+kep4GnBbRIwHbk/HneQkYBqQ3wTq\n9PpahY74Je0g3wJOqUibAkyNiKURMRuYCUwsu2KVIuK2iMjHB90HbJFed2R9k4nAzIiYHRFLgSvJ\n6tsxImJ+RDySXr8IPAGMAQ4GLk+XXQ68tz01rCZpC2B/4CJAKblj62s9czBOJE0B5kbEoxWnRgNz\nC8dzyX45O8kxwE3pdSfXdwwwp3DcSXWrImkrYBeyL7vNImJBOrUA2KxN1erJOcBngeLg3U6ur/Vg\nSE36kHQbMKqHU2eQ/Zlf7F9VD9flShkPWKO+n4uIvG/wDOC1iPhhjaw6Zfxip9SjT5LWB34MnBQR\ni6WV/x0iIiR1xGeRdCCwMCIelrR3T9d0Un2td0MqGEfEvj2lS9oBGAf8Pv3SbQE8KGkS8DQwtnD5\nFiltwPVW35yko8j+PN2nkNy2+tahsm5jWbUV3xEkrUEWiK+IiJ+k5AWSRkXEfEmbAwvbV8NVvB04\nWNL+wNrAcElX0Ln1tV64mwKIiMcjYrOIGBcR48gCxK7pz7zrgcMkrSlpHLANcH876wvZqASyP02n\nRMSSwqmOrG/yALCNpK0krUl2o/H6NtdpFcq+jS8GpkXEuYVT1wNHptdHAj+pfG87RMTnImJs+n97\nGPDLiDiCDq2v9W5ItYwbsOJPuoiYJulqsjvVy4ATojOmLZ4PrAncllrzv42IEzq4vkTEMkknArcC\nw4CLI+KJNler0h7APwOPSno4pZ0OnAlcLelYYDZwSHuq16f837pb6muJp0ObmXUAd1OYmXUAB2Mz\nsw7gYGxm1gEcjM3MOoCDsZlZB3AwNjPrAA7GQ4ikLST9NC2rOFPSuWm2WX5+alqS8yRJ20p6RNKD\nkrZusJwj06yvns7dkZbQPLCPPA5qdonNVMZuNc6fLekZSSc3k7/ZQHAwHiLSzLJrgWvTsorjgfWB\nr6bzo4C3RsROEfHfwPuAayJit4j4U4PFHUW2YFFPAvhwRNxYK4OIuCEizmqw3GIZvQ6gj4jPAt9u\nMm+zAeFgPHS8E3glIi4HSMtvfgo4RtI6ZOsjj0kLlH+BbH3cj0u6XdK6kn6WWsqPSToEQNJuqRX6\ngKRbJI2S9EHgrcAPJD0kae0e6rJi1Z30/nNTuY9JeltKP0rS+en1TyQdkV7/i6Tvp9f7Sbontd6v\nlrTeKoVIq0m6LOX7qKR/a+UP1KyVPB166NgeeLCYkFYj+wvwJuAg4MaI2AVWtKQXR8S3JH0AeDoi\nDkjnhqfujfOBgyLib5IOBb4aEcdK+gRwckQ81EtdouL1OhGxi6R3AJcAO1Zc/3+BuyXNBj4NTJK0\nCdlqe/tExCupS+PTwJcL79sFGB0RO6Z6j8CsQzkYDx19zXvvacnQPO1R4BuSziQL2Helle62B36R\n1sYYBszrI7/eTAWIiDtToF8laEbEwtRa/yXw3ohYlPqcJwD3pPLXBO6pyPePwNaSzgN+Rtb6N+tI\nDsZDxzTgg8UEScOBLcl2A+lp3WQAImKGpF2AA4CvSLoduA74Q0S8vbe39aOuy3tIewvwLKsuRn9b\nRHy4t0xS0N4JeDdwPNliOcf2o15mA8Z9xkNERNwOrFvoex0GfBO4tGIJzippZMSSiPgB8A2yP/+f\nBN4gafd0zRrK9t8DWAwMb6B6h6Y89gQWRcTiivInku2btyvwGWU7cNwH7CHpTema9SRts+rbtDEw\nLCKuBT6f3m/WkdwyHlreR7ZB6efJvoh/BnyucL6yNZsf7wicLWk52Yadx0fE0nSz7rzUrbA62fY/\n04DLgG9Lehl4e1/BHlgi6aGUxzGFsiOte/xd4KiIyIejXRIR71S2uP5USWul95wBzCi8fwxwqVZu\nyOpNOa1jeQlNK5WkXwGfiYgHC8e1bvYNVD2+SHaD8ptllmvWG3dTWNmeAy7ra9LHQJJ0NvAR4MV2\n1cGsklvGZmYdwC1jM7MO4GBsZtYBHIzNzDqAg7GZWQdwMDYz6wAOxmZmHeD/A+zW7vL/aDJ3AAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e647550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot_image(0.8, axes_units='pixels', dynamic_range=1e8, colorbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with Observed Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The observations class is compatible with the same file formats and index file that MCRE uses. You can instantiate it with a path for a directory containing the observations and the index text file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#obs = mcfost.Observations('/Users/mperrin/data/mcfost/esoha569/observations')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are ``.file_names``, ``.file_types``, and ``.file_wavelengths`` attributes that contain the entries from that table. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print obs.file_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The observations SED looks essentially identical to the model SED in its class interface. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#obs.sed.wavelength"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#obs.sed.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using these it's simple to plot an observed and model SED together, and then customize the plot however you'd like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ax = model.sed.plot()\n",
    "#obs.sed.plot(overplot=True)\n",
    "\n",
    "#ax.set_xlim(0.05, 3e5)\n",
    "#ax.set_title('Comparison of some model and data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Misc Scratch Space Code Below Here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rec.array((4000.0, 2.0, 1.0, b'lte4000-3.5.NextGen.fits.gz', True, 2.2, 0.0, 0.0, 0.0, 0.1), \n",
       "          dtype=[('temp', '<f8'), ('radius', '<f8'), ('mass', '<f8'), ('spectrum', '|S80'), ('l_is_blackbody', '|b1'), ('slope_fUV', '<f8'), ('z', '<f8'), ('y', '<f8'), ('x', '<f8'), ('fUV', '<f8')])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "par.star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(numpy.record, [('temp', '<f8'), ('radius', '<f8'), ('mass', '<f8'), ('spectrum', 'S80'), ('l_is_blackbody', '?'), ('slope_fUV', '<f8'), ('z', '<f8'), ('y', '<f8'), ('x', '<f8'), ('fUV', '<f8')])\n"
     ]
    }
   ],
   "source": [
    "names = par.star.dtype\n",
    "#values = np.asarray(par.star[0].values())\n",
    "#dtypes = ['a30' if n=='spectrum' else 'f8' for n in names]\n",
    "print(names)\n",
    "\n",
    "#np.array(par.star[0].items(), dtype=dtypes)"
   ]
  }
 ],
 "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}